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Abstract. The Lyman alpha forest power spectrum has been measured on large scales by 
the BOSS survey in SDSS-III at z ~ 2.3, has been shown to agree well with linear theory 
predictions, and has provided the first measurement of Baryon Acoustic Oscillations at this 
redshift. However, the power at small scales, affected by non-linearities, has not been well 
examined so far. We present results from a variety of hydrodynamic simulations to predict 
the redshift space non-linear power spectrum of the Lya transmission for several models, 
testing the dependence on resolution and box size. A new fitting formula is introduced to 
facilitate the comparison of our simulation results with observations and other simulations. 
The non-linear power spectrum has a generic shape determined by a transition scale from 
linear to non-linear anisotropy, and a Jeans scale below which the power drops rapidly. 
In addition, we predict the two linear bias factors of the Lycc forest and provide a better 
physical interpretation of their values and redshift evolution. The dependence of these bias 
factors and the non-linear power on the amplitude and slope of the primordial fluctuations 
power spectrum, the temperature-density relation of the intergalactic medium, and the mean 
Lya transmission, as well as the redshift evolution, is investigated and discussed in detail. 
A preliminary comparison to the observations shows that the predicted redshift distortion 
parameter is in good agreement with the recent determination of Blomqvist et ah, but the 
density bias factor is lower than observed. We make all our results publicly available in 
the form of tables of the non-linear power spectrum that is directly obtained from all our 
simulations, and parameters of our fitting formula. 
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1 Introduction 

The Lya forest is the main observational probe for studying the structure and evolution of 
the intergalactic medium (IGM). The transmitted Lya flux observed in the spectrum of a 
distant source (typically a quasar) provides us with a one-dimensional map of absorption 
along the line of sight (LOS), with the observed wavelength corresponding to the redshift of 
the intervening neutral hydrogen causing the Lya scattering (e.g. [43, 49]). 
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The way in which the flux transmission fraction is related to the underlying density 
field, temperature and peculiar velocity gradient of the gas on small scales is non-linear, and 
can be modelled in detail only from hydrodynamic cosmological simulations [12, 16, 27, 35, 
39, 44, 55, 58, 60]. However, in the limit of large scales, the transmission fraction averaged 
over a large region depends linearly on the mean overdensity and peculiar velocity gradient in 
the region, and the power spectrum of the transmission is simply proportional to the power 
spectrum of mass fluctuations, with the standard redshift distortions that were predicted 
initially for galaxy surveys [25, 29]. 

This simple linear treatment of the Lya forest applicable on large scales has led to its 
use as a cosmological tool to measure the power spectrum, first from few tens of quasar 
spectra [where only the projected one-dimensional power as a function of the parallel Fourier 
component is measured; see 16, 17, 39, 41, 58] and then in full redshift space, where the 
correlation in the transmission among parallel lines of sight is used [52]. As proposed in 
[38], and implemented in the Baryon Oscillation Spectroscopic Survey [BOSS; see 18] of the 
Sloan Digital Sky Survey-Ill [SDSS-III, see 20], the Lya power spectrum is proving to be a 
powerful tool to measure the general large-scale matter power spectrum, and in particular to 
measure the baryon acoustic oscillation scale, at a relatively high redshift which has not so 
far been probed by other observables. This measurement provides geometrical constraints on 
the expansion rate and the angular diameter distance as a function of redshift [ 8 , 19, 22, 53]. 

In contrast to the large scales, the Lya power spectrum at small scales is affected 
by a variety of non-linear physical processes governing the evolution of the IGM. These 
physical processes are highly complex, and they may include several phenomena related to 
the formation of stars and quasars in galaxies that can perturb the IGM. Some examples 
are reionization and the inhomogeneous heating caused by it, and the hydrodynamic effects 
from galactic winds and quasar jets. There is, however, a more simple assumption that can 
be made for the evolution of the IGM: that the ionization of the IGM is caused only by a 
nearly uniform radiation background, which produces a nearly uniform heating, and that 
shock waves arise only from the gravitational collapse of structure formation and not from 
the ejection of any gas from galaxies due to supernovae-driven winds or quasars. Even though 
it is known that quasar jets and galaxy winds are present in the universe and they have some 
impact on the IGM, the volume they affect may in practice be very small [40, 56, 59], and 
it is useful to test first the most simple assumption for the evolution of the IGM against the 
observations. This simple model should be mostly described by only five parameters, which 
determine the statistical properties of the Lya forest: 

• The mean transmission fraction, F(z), which depends on the intensity of the cosmic ion¬ 
izing background and is directly measured in the observations (except for uncertainties 
related to continuum fitting). 

• The density-temperature relation, usually parameterized by the two parameters To and 
7 in the power-law relation T = Td( 1 + J) 7 , where 5 = p/p — 1 is the gas overdensity. 
When the IGM is heated in photoionization equilibrium and cools adiabatically due to 
Hubble expansion, one expects this power-law relation to hold with 7 ~ 0.6, but the 
relation may be altered by the heating due to Hell reionization [28, 42, 45]. 

• The mass power spectrum of primordial perturbations near the Jeans scale of the IGM, 
A j = 2n/kj, which we can parameterize also with two parameters as a power-law with 
free amplitude and index, P(k) = A a (k/kj) na . The characteristic Jeans scale at the 
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mean density is related to the IGM temperature, although in detail it depends also on 

the entire thermal history [23], and therefore may be considered as a sixth parameter. 

Even though the large-scale properties of the Lya forest are simply understood from 
linear theory, there is a strong interest in understanding the small-scale, non-linear properties 
as well. There are several motivations for this: first, we need to test if our understanding 
of the IGM in terms of a simple uniform photoionization as mentioned above is essentially 
correct, or if there are important modifications due to a strong impact of galactic winds and 
jets [e.g., 32], or large inhomogeneities due to Hell reionization [15, 42], Second, the Lya 
forest linear power spectrum depends on two bias factors, with values that can be measured 
and can be predicted from an understanding of the small-scale physics. Finally, the Lya power 
spectrum transmission power spectrum is being measured to increasingly high accuracy, both 
in projection [from single lines of sight, see 41, 46], and in its full shape in redshift space [3], 
and a detailed comparison of the observations with predictions from numerical simulations of 
the fully non-linear power spectrum may offer us new clues to essential questions in cosmology, 
such as the impact of neutrino masses on the growth of structure [47], or limits on models of 
warm dark matter [57], or other possible variations on the nature of the dark matter. There 
is therefore a need to obtain reliable theoretical predictions for the non-linear power spectrum 
of the Lya transmission fraction as a function of redshift from numerical simulations of a 
large array of cosmological models, in terms of the most important Lya forest parameters 
mentioned above. 

The goal for the theory of the non-linear Lya forest is comparable to that of numer¬ 
ical simulations of the hot, X-ray emitting gas in clusters of galaxies. Detailed determi¬ 
nations of the gas density and temperature distributions from X-ray observations and the 
Sunyaev-Zeldovich effect, together with the mass distribution from gravitational lensing and 
the kinematic distribution of galaxies, have spurred advances in the theoretical modelling 
of clusters, the comparison of numerical codes for cosmological simulations, and tests of the 
convergence of the results. At present, the abundance of clusters of galaxies can be used to 
infer the normalization of the mass power spectrum, but this determination depends on the 
uncertain relation between the observable properties from X-rays, gravitational lensing and 
the Sunyaev-Zeldovich temperature decrement to the cluster mass. This relation needs to be 
predicted from numerical simulations, and the theoretical modelling affects the comparison 
with the power spectrum normalization derived from CMB observations [e.g., 26]. Similarly, 
the Lya forest is sensitive to the amplitude of the power spectrum and several other cosmo¬ 
logical parameters and physical properties of the IGM, but constraints on these quantities 
can only be inferred once we have a reliable understanding and modelling of non-linear effects 
on the observed properties of the Lya forest. 

The aim of this work is to study several cosmological simulations of the Lya forest 
for a variety of models, to analyze the non-linear power spectrum of the Lya transmission 
that they predict, and to test the conditions that the simulations must satisfy, in terms of 
resolution and simulation volume, to reach convergence of the results. This problem was first 
addressed in the pioneering paper of [36] (hereafter M03), and here we attempt to continue 
this study by examining a large number of hydrodynamic simulations, characterizing the 
power spectrum with a new, simpler fitting formula with several non-linear parameters, and 
studying the dependence of the linear bias factors on the IGM properties. We start in §2 by 
reviewing the definition of the linear bias factors of the Lya forest and the derivation of the 
linear power spectrum, and we introduce alternative bias factors for the Lya effective optical 
depth. The simulations and our technique for measuring and fitting the power spectrum are 
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explained in §3. The results for the non-linear power spectrum in one specific model, for 
which we have run our largest simulation, are presented in §4. Several tests of convergence 
with the box size and resolution of the simulations are performed in §5, and the results for 
the power spectrum fits for a variety of different physical models are presented in §6. Finally, 
the results are discussed in §7 and conclusions are given in §8. 


2 The Bias Factors of the Lyman Alpha Forest 


Before proceeding to describe our analysis of hydrodynamic simulations, we review here the 
standard definitions of the power spectrum and bias factors of the Lya forest. We also 
introduce a new definition for the optical depth bias factors, already discussed in the context 
of metal lines in [21], which is useful to better interpret their values and compare results at 
different redshifts and for different types of absorption systems. 

The transmission fraction F in the Lya forest is the ratio of the observed flux / that is 
transmitted from the source to the continuum flux, f c , when there is no intervening absorp¬ 
tion. Usually a model of the source continuum spectrum is used to calculate F = f /f c from 
the observed flux. The fluctuation in the transmission is defined as 


Sf 


F 

W) 


( 2 . 1 ) 


where F(z) is the mean value of the transmission, evaluated as a function of redshift. The 
power spectrum Pp(k, )i) is that of 5f in redshift space, with /i being the cosine of the angle 
of the Fourier wavevector of modulus k from the LOS to the observed source. On small 
scales, the value of 5f at a certain point depends on the complex non-linear evolution of the 
IGM. However, the average of Sf over a large enough scale is related to the mass density 
fluctuations according to a linear expression of the form 5p = bp$6 + bp r] r). Here, 5 is the 
mass density fluctuation, and 


V = 


1 dv p 
aH dx r , 


( 2 . 2 ) 


is the dimensionless gradient of the peculiar velocity v p along the LOS, both averaged over 
the same large region that is used to average 6f ; x p is the comoving LOS coordinate, 77 
is the Hubble constant, and a the scale factor. This is the most general linear expression 
for any tracer of large-scale structure. The reason is that the only scalar quantities that 
can be constructed from the deformation tensor, (j)jj (equal to the second derivatives of the 
gravitational potential), and the LOS unit vector, m, are the trace of </>, which is proportional 
to 6, and which is proportional to r\. The bias factors are simply defined as the 

partial derivatives of bp with respect to <5 and rj, after the large-scale averaging is done: 


bps = 


dS F 


bFn = 


d6 F 

dr) 


(2.3) 


where each partial derivative is understood to be done by holding the other variable (5 or rj) 
constant. 

The bias factors in equation (2.3) are not defined in the same way as for point objects like 
galaxies. In general, denser regions or regions where the Hubble expansion rate has slowed 
down have stronger absorption, and therefore lower transmission F, so the values of bps and 
bp p are negative. Moreover, when F is close to unity (which occurs at low redshift), the 
absolute values of the bias factors are very small simply because they express the absorption 
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fluctuation compared to the total transmitted fraction, when the mean absorption is very 
small. A physical bias factor should reflect the fluctuation of a quantity that is zero when the 
mass density is zero, whereas the transmitted fraction is obviously one when the gas density is 
zero. The quantity that reflects the relative amount by which the Lya absorption fluctuates 
when the mass or the peculiar velocity gradient fluctuate is obtained from the fluctuation of 
an effective optical depth, r e = — log F (where F has been averaged over a large, linear scale 
before taking the logarithm), which has a relative fluctuation 5 r = 5 f/ log F = —5p/f e . The 
bias factors for this effective optical depth fluctuation are: 


_ d5 T _ bps 
85 log F 


85 T b Fr/ 

Tr> dr] log-F 


(2.4) 


The usefulness of these new definitions is made apparent by considering various simple 
models for the Lya forest. Let us first imagine that the Lya absorption systems are clouds 
of gas distributed in space with a density bias factor b c \ in other words, wherever the mass 
density fluctuates by 5, the number density of gas clouds fluctuates by b c 5, but their individual 
absorption line profiles do not vary. If the internal dynamics of these clouds are not aligned 
in a correlated way with the LOS depending on the value of ?y, the effective optical depth 
they produce in the spectrum fluctuates as 5 r = b c 5 + r), because a peculiar velocity gradient 
is simply squeezing the absorption lines of the clouds in the spectrum by the factor 1 — 77 , 
without altering the absorption line profiles. Hence, for these population of clouds we have 
b T s = b c and b Tr) = 1 : the density bias factor reflects the true physical bias of the population 
of clouds, and the peculiar velocity gradient bias factor is equal to one, just like for any 
population of objects that is selected isotropically, i.e., independently of the LOS direction. 

Next, consider the case where the IGM absorption is optically thin everywhere. The 
optical depth is proportional to the density of hydrogen atoms, which on large scales must 
again behave like a biased population of objects with some bias factor b a , and b TV must 
again be unity because the column density of hydrogen atoms in a given interval of the Lya 
spectrum is proportional to the integrated optical depth, which cannot be modified by any 
shifts due to peculiar velocities. 

Therefore, if we consider that the true Lya forest is a combination of optically thin 
absorption by the IGM, plus a population of clouds with optically thick absorption lines that 
follow a linear bias factor but do not change their internal properties as a function of 6 and 
do not have correlated orientations with the principal axes of the large-scale deformation 
tensor, we conclude that b T . q = 1 , and b T s reflects a true, physical bias factor that results 
from a weighted average of the optically thin, intergalactic neutral hydrogen bias b a , and the 
gas clouds bias b c . 

This model is of course not exactly correct, because there is an intermediate density 
range of gas that is not optically thin and is not in a population of clouds that have lost any 
alignment of their internal dynamics with the surrounding large-scale structure. Nevertheless 
we can expect it to provide a first approximation to the reality of the Lya forest, and in this 
way the value of the two bias factors can have a physical interpretation. Note also that these 
optical depth bias factors can be defined for any other set of Lyman series lines and metal 
lines, simply by using the appropriate value of F in equation (3.1), as proposed in [21], and 
then the bias factors of any hydrogen or metal absorption lines associated with a certain 
population of clouds or galaxies are equal to the true bias factors of these objects. 

The linear power spectrum is derived as in [29], from the simple fact that in the linear 
regime, r] = f(Q)fi 2 5 in Fourier space, where f(Fl) = dlogG/dloga, and G{a) is the growth 
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factor. For the linear power spectrum of the transmission fluctuation 5f, we have 


Pf( k, n) — + fill 2 ) 2 Pp{k) + N 0 ) 


(2.5) 


where the redshift distortion parameter is 


bFt] f (Q) 
bFS 


( 2 . 6 ) 


and Pp is the matter density fluctuation power spectrum. We have included the intrinsic 
shot noise term JVo, which is a shot noise that should remain even when Pp is measured from 
an arbitrarily dense set of absorption sightlines due to the random nature of the formation 
of absorption line systems as the evolution of the IGM becomes non-linear [37], although 
this is believed to be very small and we shall ignore it in this paper. This linear power 
spectrum is of course valid only in the limit of large scales, and the models we shall use to 
fit our simulations include a non-linear multiplicative term, as described below in equation 
(3.3. The linear power spectrum of the effective optical depth 5 r is the same, except for 
the normalization, which changes by replacing bps by b T $. Note that (5 is the same for the 
transmission or effective optical depth power spectra, because the ratio of the two bias factors 
remains unaltered. 

Finally, we mention also the radiation bias factor bp p, defined as the variation of 5f 
when the photoionization rate of hydrogen T, determined by the intensity of the ionizing 
background, varies. If dr is the relative fluctuation of this photoionization rate, then the 
total transmission fluctuation averaged on a large scale is 5f = bps5 + bp v r] + bppdr. The 
modification of the power spectrum Pp due to fluctuations in the ionizing radiation intensity 
caused by sources that are tracers of the mass density fluctuations was discussed in [24, 48]. 
We ignore here the possible additional effect from Hell reionization of a large-scale modifi¬ 
cation of the temperature-density relation. Under the assumption that the intensity of the 
ionizing background does not appreciably affect the temperature and hydrodynamic evolu¬ 
tion of the IGM, and changes the optical depth at every spectral pixel in inverse proportion 
to the photoionization rate, then the variation of the transmission fluctuation with dr can be 
computed in terms of the probability distribution of the transmission fraction P(F), as 


5f 


i [ dF P(F) exp 

L Jo 


f log F \ 

\l + 5r) 


Jo dF P{ F) F log F 
- p - d r 


feFrdr • 


(2.7) 


Just as before, we define the optical depth radiation bias factor as 


b FT _ Jo dF P(F) F log F 
log F F log F 


( 2 . 8 ) 


An alternative simple model for computing the peculiar velocity gradient bias factor is 
to assume that all the fluctuations determining the Lya forest absorption spectrum can be 
treated in the linear regime. Using this assumption, [51] showed that bp^ should be given 
by the same expression for bpp in equation (2.8). The reason is easy to understand: for 
linear fluctuations, the optical depth at any spectral pixel is simply multiplied by the factor 
(1 — r ?) -1 ~ I + 77 under the effect of a peculiar velocity gradient, so the same derivation shows 
that bp v = bpp. For the case of radiation fluctuations, however, the derivation of equation 
(2.8) does not need to assume that the Lya forest fluctuations are linear even on small scales, 
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so the prediction for the radiation bias factor is much more reliable. We will show in §4 that 
in fact, bp v is quite different from bp r because of the non-linearities that affect the change 
of small scale Lya forest fluctuations under a variation of the large-scale peculiar velocity 
gradient. 

3 Method of Analysis of the Simulations 

Our goal in this paper is to use cosmological hydrodynamic simulations of the IGM to predict 
the three-dimensional power spectrum of the Lyman alpha forest in redshift space, Pp(A;, /j; z), 
where k and /x are the modulus and the cosine of the angle from the LOS of the Fourier mode 
vector, and z is the redshift. The simulations used are described in §3.1. The method of 
analysis is inspired in that of M03 and is based on the following steps: (1) starting from a grid 
of cells containing the hydrodynamic quantities of gas density, ionized fraction, temperature 
and velocity at a certain redshift output of a simulation, the corresponding spectra of Lya 
transmission are computed for the entire grid, using one of the simulation axes as the assumed 
LOS direction, and the three-dimensional Fast Fourier Transform of this transmission field is 
obtained, as described in §3.2; (2) the mean value of Pp(L,/x) is computed in bins of (fc,/x), 
and errorbars are assigned which take into account the variance due to the finite simulation 
volume (§3.3); (3) a parameterized fitting function for Pp(k,n) is chosen to obtain best-fit 
values of the parameters for several simulations (§3.4). 

3.1 Simulation characteristics 

Two types of hydrodynamic simulations will be used in this chapter. Most of the simulations 
rely on the Tree-PM (Particle Mesh) Smoothed Particle Hydrodynamics (SPH) GADGET- 
II code [54], and the bulk of our analysis will be performed on the outputs of these Lagrangian 
simulations. One simulation that is based on a fixed-grid Eulerian code is also used, to allow 
for a first comparison of the results for the two types of hydrodynamic numerical methods. 
This simulation is described in [11], 

Table 3.1 shows a list of all the simulations that will be used in this work, including 
variations in the spatial grid size and the spectral pixel size for the analysis of the Lya forest. 
The first two columns give the comoving box size, L, and the number of dark matter particles 
in the SPH simulations (the number of dark matter and gas particles in the SPH simulations 
being the same). The third one gives the number of cells, N c , in the uniform grid that is 
constructed to compute the density, temperature and velocity in real space. Note that the 
simulation labelled Euler does not use particles, and the cells used to run the simulation are 
N% and are directly used as this spatial grid. The Lya spectra are computed for each of the 
three axes of the simulation playing the role of the LOS, with the number of pixels in each 
spectrum from each row of N c cells of length L given in the fourth column; generally there 
are as many pixels in the spectra as cells in the spatial grid, except in the analysis labelled 
P1024 where the number of pixels is doubled. The other columns give values of physical 
parameters used in the simulations: the variable as parameterizing the present amplitude of 
linear perturbations on a sphere of 8 /i _1 Mpc, the mean temperature at the mean density To, 
and the power-law index that fits the density-temperature relation at low densities, which is 
described below in more detail. In general, models have variations of different parameters 
around the values of the fiducial model in the first row of Table 3.1, and they are labelled 
with names that refer to the parameter that is being varied. 
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Name 

Box size 

Particles 

N c 

Pixels 

08 

7 

log (To) 

Fiducial 

60 Mpc/h 

512 3 

512 3 

512 

0.8778 

1.6 

4.3 

P1024 

60 Mpc/h 

512 3 

512 3 

1024 

0.8778 

1.6 

4.3 

C256 

60 Mpc/h 

256 3 

256 3 

256 

0.8778 

1.6 

4.3 

R384 

60 Mpc/h 

384 3 

512 3 

512 

0.8778 

1.6 

4.3 

R384C 

60 Mpc/h 

384 3 

256 3 

256 

0.8778 

1.6 

4.3 

R640 

60 Mpc/h 

640 3 

512 3 

512 

0.8778 

1.6 

4.3 

R640C 

60 Mpc/h 

Oi 

O 

CO 

256 3 

256 

0.8778 

1.6 

4.3 

L80 

80 Mpc/h 

512 3 

512 3 

512 

0.8778 

1.6 

4.3 

L120 

120 Mpc/h 

768 3 

512 3 

512 

0.8778 

1.6 

4.3 

Euler 

50 Mpc/h 

— 

2048 3 

2048 

0.82 

1.54 

4.03 

Lagrange 

50 Mpc/h 

512 3 

512 3 

512 

0.82 

1.58 

4.10 

Planck 

60 Mpc/h 

512 3 

512 3 

512 

0.8338 

1.4 

4.2 

G1.3 

60 Mpc/h 

512 3 

512 3 

512 

0.8778 

1.3 

4.3 

G1.0 

60 Mpc/h 

512 3 

512 3 

512 

0.8778 

1.0 

4.3 

G1T4 

60 Mpc/h 

512 3 

512 3 

512 

0.8778 

1.0 

4.0 

SO. 76 

60 Mpc/h 

512 3 

512 3 

512 

0.7581 

1.6 

4.3 

SO.64 

60 Mpc/h 

512 3 

512 3 

512 

0.6396 

1.6 

4.3 


Table 1 . List of simulations and analysis variations used in this paper. The first column lists the name 
we give to the simulation/analysis, the second the simulation box size, and the third indicates the 
number of particles used in the simulation (both for dark matter and gas, except for the simulation 
named Euler, which uses a fixed Eulerian grid instead of gas particles). The fourth column gives 
the number of cells, N c , used to represent the hydrodynamic variables in the spatial grid that is 
computed to obtain the Lya forest spectra, and the fifth column is the number of pixels on the line 
of sight direction used to compute the Lyct spectra. The last three columns give the power spectrum 
amplitude and the (To, 7 ) parameters of the temperature-density relation in the simulation, where To 
is expressed in kelvin. 


3.1.1 SPH simulations 

All simulations in Table 3.1 except for the one denoted as Euler were run using the publicly 
available Tree-Particle Mesh Smoothed Particle Hydrodynamics (SPH) GADGET-II code 
[54] ' 

The fiducial simulation uses a box of 60 comoving h 1 Mpc and 2 x 512 3 particles (for 
the total of gas and dark matter). Other simulations are run with larger boxes of 80 and 120 
h _1 Mpc (L80 and L120) to test the effect of the missing large-scale power, or with different 
resolution to check the convergence as the particle masses are reduced. The cosmological 
model is flat ACDM with the following parameters, using standard notation: Ho m = 0.3, 
A = 0.7, Hob = 0.05, Hq = 70kms~ 1 Mpc _1 , n s = 1 and erg = 0.8778. The initial conditions 
are generated using the software CAMB 1 and the Zel’dovich approximation at the initial 
redshift of z = 49. The particle mesh grid used to calculate the long range forces is chosen to 
have the same number of cells as the number of gas particles, 512 3 , while the gravitational 
softening is 4 kpc /h in comoving units for the 60 h -1 Mpc box and scales proportionally to 
the initial particle separation for the other simulations. The hydrodynamical processes are 
followed according to the prescription of [30]. Star formation is also included in the model 
with a simplified prescription that allows to instantaneously convert into a star particle any 

1 http: //camb .info / readme .html 
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gas particle of overdensity larger than 1000 and temperature colder than 10 5 K. This has 
been demonstrated to have negligible impact on the Lya forest transmission power spectrum 
[58] ‘ 

Most of the SPH simulations in Table 3.1 are based on the same cosmological model as 
the fiducial one. The exceptions are the simulations SO.76 and SO.64, where the amplitude of 
the initial power spectrum is varied, and the Planck and Lagrange simulations. The Planck 
simulation uses a model that is consistent with the most recent CMB measurements from the 
Planck mission (REFERENCE), with the following parameters: Ho m = 0.3175, A = 0.6825, 
Hob = 0.049, Ho = 67.11 km s _1 Mpc -1 , n s = 0.9624 and erg = 0.8338. The Lagrange 
simulation is run for the same cosmological model as the Euler simulation that is run with 
the Eulerian code, described below. 

Our main goal in analyzing these simulations is to understand how the non-linear 
power spectrum varies several physical parameters, such as the amplitude erg or the density- 
temperature relation. These comparison are limited by the intrinsic random variations in 
the measured power due to the sample variance in simulations of limited box size. To reduce 
this intrinsic sample variance, we have chosen to generate the initial conditions of all the 
GADGET-II simulations by setting the amplitude of every Fourier mode to the exact rms 
amplitude predicted by the power spectrum, instead of generating it following the Rayleigh 
distribution. The mode phases are still generated randomly. Only for the fiducial simula¬ 
tion, we have generated several realizations with different random seeds, including some cases 
where the Rayleigh distribution for the amplitudes is included, which will be analyzed in §5.2 
to check for any effect that this can have on our results. This implies that the power spectrum 
of the initial conditions for our simulations is exactly equal to the value predicted by the cos¬ 
mological model at each Fourier mode, without any variations due to sample variance. The 
Lya transmission power spectrum that is obtained, however, has random variations caused 
by non-linear Fourier mode couplings. 

The impact of different thermal histories on the Lya forest is explored by modifying the 
Ultra Violet (UV) background photo-heating rate in the simulations, as in [ 6 ]. A power-law 
temperature-density relation, T = To(l + 5) 7_1 , arises in the low density IGM as a natural 
consequence of the interplay between photo-heating and adiabatic cooling [28]. Two different 
values for the temperature at mean density, To, and three different values for the power-law 
index of the temperature-density relation, 7 , are considered to examine the impact of the 
temperature-density relation on the Lya power spectrum; the most recent observational 
constraints [4] favor values for these parameters close to those of our fiducial model. The 
different thermal histories are constructed by modifying the fiducial simulation He II photo¬ 
heating rate according to tueii = a x e'j id HeII , changing the parameters a and u [ 6 ]. The 
fiducial thermal history appears to be in overall good agreement with recent determinations 
based on line profile fitting [5]. The thermal histories for the SPH runs are built to guarantee 
that the p — T relation is approximately constant with redshift in the range z = 2.2 — 3, so 
the quoted values of 7 and To approximate reasonably well the density-temperature relation 
in the whole redshift range investigated here. 

3.1.2 Eulerian simulation 

In addition to the SPH simulations based on a Lagrangian approach, a simulation based on 
an Eulerian code is used, [described in 9, 10, 13, 14], The simulation used here is on a 50 
Mpc/h box with 2048 3 cells, and uses the cosmological model with parameters Ho m = 0.28, 
Hob = 0.04, A = 0.72, Ho = 70kms _1 Mpc -1 , n s = 0.96, and <7g = 0.82. This simulation has 
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been run with standard initial conditions, with the amplitudes of every Fourier mode following 
the Rayleigh distribution. The p — T relation is determined by following the photoionization 
heating that is derived from a model of the ionizing background that is computed as the 
simulation is run. In the density range that is important for the Lya forest over the redshift 
range 2.2 < z < 3 that we analyze in this paper, this p — T relation is well approximated by 
the values of To and 7 in Table 3.1. 

To compare this simulation with an equivalent one run with the SPH method, we have 
run the simulation designated ’’Lagrange” in Table 3.1. The Lagrange simulation is run for 
the same cosmological model and box size as the Euler one, and using Hell heating parameters 
to approximately mimick the p -T relation in the Euler one. This is further discussed in §5.3. 

3.2 Extracting the Lya power spectrum from the simulations 

We now start discussing the full procedure for processing the simulation outputs to obtain 
an estimate of the Lya transmission power spectrum, and for analyzing fits to this power 
spectrum. This procedure is summarized in the diagram in figure 1, and discussed in detail 
in the rest of this section. 

For the GADGET-II simulations, the SPH formalism for computing the hydrodynamic 
variables of gas density, temperature and velocity on a Cartesian grid, and then extracting 
mock Lya spectra, is followed as described in the Appendix A4 of [55]. For the Euler 
simulation, the Cartesian grid that the simulation is run on is used directly to obtain the 
Lya spectra, as in [44]. Each simulation grid in real space is used to generate three distinct 
boxes of Lya spectra, taking each of the three axes as the LOS. For each of the three axes, 
the spectra for the entire simulated box are computed, resulting in Lya spectra. 

Apart from the parameters of each simulation, an additional parameter is necessary to 
compute the Lya spectra: the intensity of the ionizing background, which can be altered to 
adjust the mean transmission F(z) to a certain value. The mean transmission fraction is 
fixed to the value given by the expression 

F(z) = exp [-0.0023(1 + z) 3 ' 65 ] , (3.1) 

which was found to adequately fit the observational data of high-resolution spectra by [31], 
after subtracting the estimated metal contribution. We note that more recent determinations 
give comparable values of F. but there are substantial uncertainties in this determination 
[1,4]. 

The computed Lya spectra are modified to adjust this value of the mean transmission 
by using the approximation that the optical depth varies at each pixel as the inverse of the 
intensity of the ionizing background, and that the gas temperature is not affected by this 
background intensity. This assumes that collisional ionization can be neglected and that the 
atomic fraction is much smaller than unity, which is generally an excellent approximation 
(except in high density regions where the optical depth is very large in any case, and therefore 
does not affect the computed Lya spectra). The use of this approximation avoids having to 
recompute the Lya spectra every time that the mean transmission is adjusted to the required 
value in the expression above. The assumption that the temperature and hydrodynamic 
evolution of the IGM is not affected by the intensity of the ionizing background is not as 
accurate because cooling by line excitation is neglected, but in any case, here we are interested 
in examining the dependence of the predictions for Pp on an assumed, fixed p -T relation, 
and separately on F. 
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Figure 1. Flow chart of the analysis method of the power spectrum from simulations followed in 
this work. A set of GADGET-II simulations are used with different resolution and box sizes, and 
different physical properties: mean transmission F, power spectrum amplitude as, mean temperature 
To, and temperature-density relation slope 7 . The output of the simulations at various redshifts is 
represented in a real space grid (the Eulerian simulation output directly provides this grid), and then 
the Lya optical depth is computed in redshift space. These optical depths are rescaled to keep the 
mean transmission fraction fixed, and then the Fast Fourier Transform is performed and the square 
moduli provide estimates for the power spectrum Pp(k x , k y ,k z ). The power spectrum is then averaged 
in bins of (k, p), errors are estimated, and fits to our proposed equation (3.6) are done. At first, tests 
of convergence with parameters related to the binning of Fourier modes and the errors are carried out. 
After these tests are performed, a fiducial model is chosen and further convergence tests are done for 
the simulation box size, resolution, and grid cells. Finally, the dependence of the results on physical 
effects is examined. 


We will generally present results at the redshifts z = 2.2, 2.4, 2.6, 2.8, and 3, where 
equation 3.1 implies rescaling the Lya spectra to mean transmission values of F = 0.8517, 
0.8185, 0.7813, 0.7404, and 0.6960, respectively. We shall use these values except for a few 
cases discussed in section 6.3, where we examine the variation of Pp under changes in F. 

A Fast Fourier Transform is applied to the entire box of Lya spectra, for each of the 
three cases taking each axis as the LOS. Usually, the grid is cubic with N% cells, except in 
the P1024 model in Table 3.1 where the number of pixels is 2 N c . For the latter case, we 
first average the value of F = exp(—r) in every two pixels to obtain a cubic grid, and then 
compute the Fourier transform. The routine fft from the scipy package in Python is used 
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for this computation 2 . This results in Nj}/2 independent Fourier modes for each of the three 
axes chosen as the LOS, each one with a modulus and a phase. The moduli are used to 
obtain the estimate of the Ly a power spectrum Pp(fe, //; z). 

3.3 Power spectrum estimation, Fourier space binning, and errorbars 

We discuss in this subsection the procedure for evaluating the power spectrum in Fourier 
space bins, using the Fourier modes of the Lya transmission field of the simulations. The 
wavenumbers that are available in a simulation of box size L go from the minimum value, k\ = 
27 t/L, to the maximum value equal to the Nyquist frequency, k\N c /2. Several parameters 
are involved in our choice of binning in Fourier space and the computation of error bars of 
Pp used to obtain fits, which we enumerate here: 

1. The cutoff scale k c . This is the maximum scale at which we consider that Pp can be 
reliably predicted from a simulation and measured from the observations. We use only 
modes with k < k c to fit the results of Pp, where k c is less than the Nyquist value. For 
our fiducial simulation, the comoving cell size is L/N c = 0.12 h~ l Mpc, corresponding 
to a velocity width ~ 12kms _1 and a Nyquist value ttN c (1 + z)/(HL) ~ 0.25skm - . 
Most of the modes in a simulation have wavenumbers near the Nyquist value, which 
are affected by the absorption tails of high column density systems arising from highly 
non-linear collapsed structures that may not be correctly modelled in the simulations. 
Moreover, in practice the observed power at this large k is highly sensitive to the 
presence of narrow metal lines, which are not included in the simulations. Following 
M03, we consider that any comparison of theoretical and observed power spectra is not 
reliable for wavenumbers above k( 1 + z)/H ~ 0.1 s/krn, which is approximately a fixed 
comoving scale in our examined redshift range We therefore choose, for our fiducial 
simulation, k c = 100 k\ = 10.47 L/Mpc, and for all simulations we keep the physical 
value k c = 10.47 h /Mpc fixed. All the modes with k > k c are discarded. This leaves, 
for the fiducial model, 4 x 10 6 independent Fourier modes with k < k c to be used in 
our analysis. 

2. The transition scale kt . For effectively computing a y 2 function to fit the estimated 
power spectrum from a simulation to an analytic model, having more than 10 6 values 
from independent Fourier modes is still a very large number, and the vast majority of 
these Fourier modes are at high k. To reduce this number, we define bins in the ( k , /r) 
variables to average the estimated power spectrum within each bin when k is larger 
than a transition scale kt- For k < kt, different modes obtained from the simulation are 
averaged only when they have exactly the same values of ( k,n ), and compared to the 
power spectrum values from an analytic model at exactly the same ( k, fi) to compute 
the x 2 function. 

We choose kt = 1 L/Mpc for the fiducial simulation with box size L = 60 h -1 Mpc, 
which is close to the geometric average of k c and k\. This results in roughly the same 
number of bins at k > kt, as different values of ( k,[i ) at k < kt, which optimizes the 
efficiency and accuracy of the calculation. For different box sizes, we keep fixed the value 
of ktL = 60, so that the number of Fourier modes evaluated without binning remains 
roughly constant. For example, for our largest simulation with L = 120 h -1 Mpc, we 
use kt = 0.5 h/ Mpc. 

2 http: //docs.scipy.org/doc / numpy / reference / routines.fft .html 


12 - 



3. Number of bins in (k,n). The range of k from kt to k c is divided into 16 bins that 
are equally spaced in log fc, and /x is divided also into 16 linearly spaced, equal bins 
from 0 to 1. This gives a total of 256 bins for k > kt for which the power spectrum is 
estimated, for all simulations. Simulations with larger box size, with a smaller value of 
kt , have therefore a larger bin size in log k. The number of different values of (k. /x) that 
are obtained at k < kt turns out to be 296 for our choice of ktL = 60, therefore a total 
of 552 fitting points are used for all simulations when fitting Pp(k,n) to an analytic 
model. 


4. 


Errors ap(k,fi). Evaluating the x 2 function requires assigning an error to each Fourier 
mode estimated from a simulation. We assume that the Fourier modes are independent, 
and therefore that our covariance matrix is diagonal. We keep a count of the number 
of Fourier modes, n F (k,fi), contributing to each of the 552 values of (/c,/x) (including 
all three projections we use of the simulation box in redshift space, so n F is always a 
multiple of 3). We then compute the error of each power spectrum evaluation according 
to: 


vp{k,n) 


P F {k,fi) 


1 /y/n F (k,fj,) + e 


(3.2) 


where Pp is the estimated value of the transmission power spectrum. For e = 0, this 
is the expected error owing to the Poisson variance associated with the number of 
independent modes available in the simulation. The constant e is included to avoid an 
excessive weight to the overall fit from the modes with the highest values of k, following 
the procedure of M03. The number of modes below a certain value of k grows as k 3 , 
so if a reasonable value of e is not included, any analytic fit will need to be very highly 
accurate for all the high-k modes before the low-k modes are of any importance in 
determining the minimum of the x 2 function. Unfortunately, there is no clear objective 
way to decide the value of e that one should choose to obtain a fit that is adequately 
weighting the results of a simulation over the broad range of k that is being probed, 
and the results of the fits depend on e. In our work we have chosen a constant value 
of e = 0.05 for all the simulations and analyses. The optimal value of e is discussed in 
Appendix A. 


It is worth going through some examples of the number of independent modes available 
for the power spectrum estimate for the smallest wavenumbers, starting at k\ = ( 2ir)/L . For 
the direction parallel to the LOS (with k x = k y = 0 and k z = k\ , choosing the z-axis as the 
LOS), with /x = 1, only one independent Fourier mode is obtained (the mode with k z = —k\ 
is not independent because of the condition that the Lya transmission held is a real function). 
Two independent modes perpendicular to the LOS, with /x = 0, are available, for k x = k±, 
k y = 0, and k x = 0, k y = k\. The next smallest modes have k = y/2k\, with two independent 
modes for /x = 0 (k x = k\, k y = Aq, and k x = k \, k y = — k\), and four independent modes 
for fi = l/y/2 (with either k x or k y being equal to +/- k z ). In general, eight independent 
modes are available for any values of k x , k y and k z when they are all different from zero and 
k x 7 ^ k y (owing to the symmetry under two independent sign changes and under the exchange 
of k x for k y ), which are used to estimate the Lya power for the values k = (k 2 + k 2 + k 2 ) 1 / 2 , 
/x = k z /k. In this case, the estimate of the power at this (/c,/x) will come from n F = 24 
values, because eight independent modes are obtained for each of the three axes chosen as 
the LOS. For some modes, this number is further increased whenever several combinations 
of (k x ,k y ,k z ) yield the same values of (/c,/x) (e.g., for k x = 3, k y = 4 and k x = 5, k y = 0). 
There are 296 different values of (k, /jl) that are obtained from a cubic box with kL < 60. 
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For the bins at k > kt, the average values of (k. /x) of the contributing Fourier modes in 
each bin are stored, in addition to the mean power spectrum value. These average values are 
usually very close to the central values of the bin (because the number of modes included in 
each bin is large), but they are not exactly equal. The model of the power spectrum to be 
fitted is then evaluated at these average values instead of the bin center. 

3.4 Parameterized fitting function for the Lya power spectrum 

The Lya transmission power spectrum obtained from the simulations will be fitted to the 
following analytic model: 


P F (k, /x) = b 2 F5 (1 + /3/tx 2 ) 2 P L {k) D(k, /i) . (3.3) 


The first terms on the right hand side are derived from the linear perturbation theory of [29] , 
as explained in §2, where PrXk) is the mass density fluctuation linear power spectrum, bps 
is the density bias factor of the Lya transmission, and /3 the redshift distortion parameter. 
The function D(k,fd) is the deviation from linear theory due to non-linear evolution, so we 
expect D to approach unity in the limit of small k. 

We shall use two different fitting models for D(k,/i ) in this paper, although we have 
tested many others before deciding on a formula that provides good fits. First, the expression 
used by M03, which we designate Dq, with a total of 8 free parameters, 


D 0 (k,n) = exp 




k ■ fj, 


k v o(l + k/k v i) a " 1 


CkvO~ 


(3.4) 


Second, the expression we shall use in most of our fits, D i, is a new one that has only 6 
free parameters, and actually only 5 are used in most of our fits. We make the ansatz that 
the non-linear correction should behave as D — 1 oc k^PL(k) in the limit of small k, because 
perturbation theory predicts that the second order terms of Pp should depend on integrals of 
products of four linear perturbations. First, we define the Fourier amplitude of linear density 
fluctuations as 

A \k) = ±k 3 P L (k ) . (3.5) 

The fitting formula D\ that we use is 


D\(k, n) = exp 


[qiA 2 (k) + q 2 A i {k )] 


1 




(3.6) 


These equations are to be understood as simple fitting formulae that have been found to pro¬ 
vide useful fits to the numerically obtained power spectra from the simulations by experience. 
However, some physical motivation for the various terms can be provided as follows: 

• Non-linear enhancement: The power spectrum is increased on scales near the onset 
of non-linearity, because non-linear collapse of structure tends to enhance the power 
relative to the linear prediction. For the D\ formula, and ignoring for now the /x de¬ 
pendence, this term is q\A 2 [k ), with the optional addition of the higher order term 
q 2 A‘ l (k) that may be included to improve the fit. The constants q\ and q 2 are dimen¬ 
sionless and control the importance of this non-linear enhancement. In the Dq formula 
this term is ( k/k n i ) anl . The scale k n i at which non-linear effects start being important 
should roughly obey A 2 {k n i) ~ 1, so the dimensionless constants q\ and q 2 in the D\ 
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formula are expected to be of order unity. The local slope of the power spectrum near 
this non-linear scale at z — 2.5 is n e ff — —2.3, implying that A 2 ~ k o.7 

near this scale, 

which agrees with the typical value found for a n i in M03. In the limit of small k , the 
difference D — 1 is proportional to A 2 (k) in our new formula, so an extrapolation to k 
values smaller than those probed by our simulations can give a reasonable prediction, 
while the formula Dq of M03 is expected to overpredict D — 1 for small k. 

• Jeans smoothing: The gas pressure suppresses the power below the Jeans scale. Small- 
scale power is present in regions of high density, where the Jeans scale is reduced in 
the highly non-linear regime, but this should not greatly affect the Lya forest which is 
mostly sensitive to moderately overdense structures. This power reduction is modelled 
by the isotropic term with the scale k p playing a role that is reminiscent of a Jeans 
scale, although this cannot be taken literally because our formulae are simply a fit to 
a highly non-linear numerical result. The Dq formula has a free power-law a p for this 
term, but we have found that good fits are obtained by fixing a p = 2, and therefore we 
fix this in the D\ formula. 

• Line-of-sight broadening: Finally, non-linear peculiar velocities and thermal broadening 
cause a smoothing of the correlation along the LOS, and therefore a suppression of 
power that increases with ^ i. This suppression was found in M03 to be well matched by 
a power-law dependence on n inside the exponential, but several parameters had to be 
added to fit the ^’-dependence of this non-linear anisotropic term, the third one in the 
Dq formula. We find that by multiplying this term by the same non-linear enhancement 
depending on A 2 (k), a simple power-law dependence on k also provides a good fit. 

Therefore, our new formula D\ has two advantages over Dq: it has the correct behavior 
for D\ — 1 in the limit of small k , and it reduces the number of non-linear parameters from 
eight to six, and in fact to five in most of the simulations analyzed in this paper where we 
will set q 2 = 0, still providing sufficiently good fits. 

To calculate the fit of the values of Pp(k, /r) extracted from the simulations to equation 
(3.3), we use the Montecarlo Markov Chain method [hereafter, MCMC; see, e.g., 2, for a 
review on the Metropolis algorithm that we use] to minimize the % 2 function, which we 
compute as 

2 _ (Ps — Pm ) 2 _sr^ (Ps/Pm —Ps) 2 / ? ? \ 

x 2^ {apPm /p s) 2 Z. a 2 p ’ v-n 

where P s is the power spectrum measured from the simulation, P m is the power spectrum 
of the model being fitted, and the sums are over all the bins in (fc,/i). The errors ap are 
computed using equation (3.2) with Pp = P s , so that they do not depend on the model, 
and are the ones shown in our figures. However, the % 2 function is obtained with the errors 
computed from the fitted model, which are ap P m /P s . 

4 Results: the L120 Simulation 

This section presents the results of the non-linear Lya transmission power spectrum for the 
L120 simulation, the largest SPH simulation we have analyzed, with 768 3 particles (see Table 
3.1). We start using the fitting formula D i in equation (3.6), which will be used in all our 
results for other simulations. Figure 2 shows the estimated ratio Pp(k, fi)/Pp(k) at z = 2.2 
as colored points with error bars, and the resulting fits as curves, in the left panel. The 
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non-linear term D(k,/i) is shown in the right panel. The colored curves are the fit to the 6 
free parameters in D\{k ,//) plus the two linear bias parameters, and the black curves show 
the fit that is obtained when the q 2 parameter, multiplying the second-order term in A 2 (A:), 
is set to zero. 

The structure of these two panels will be the same for the various models analyzed in 
the rest of this paper. A set of four curves and points are shown for each fit, corresponding 
to the four intervals of /i indicated in the figure. An additional two curves show in this case 
the fit model at /i = 0 and /x = 1. The points are obtained by averaging the modes within 
the bins in k that are plotted and these four bins in /x. Even though the fit is performed 
with 16 bins in /x at k > kt , as described in §3.3, the results are then further averaged into 
4 bins for the purpose of display only. The fit also uses, as described above, 16 bins in log k 
within kt < k < k c , which are directly plotted, and individual mode values for k < kt- These 
individual mode values are also averaged, for display purposes, into 16 bins in log k between 
k\ = 27 t/L and kt , allowing for easy visualization of the results in plots that are similar to 
those in M03. Note that, for small k, some of the bins in log k and /x do not include any of 
the actual values of (fc, /x) from the simulation modes, and in this case they are absent from 
the plot. 

The way in which the values and errors of the power spectrum estimates in the original 
bins used for the fit are averaged into the bins used to make the figures is as follows: each 
original bin is assigned a weight w{k, /x) = l/a 2 P (k, /x), using the errors in equation (3.2). The 
values of the power spectrum and the mean coordinates of the coarse bins for the plot are 
obtained by averaging Pf/Pl, log k and /x with these weights, and the new error in the coarse 
bin is set to the inverse square root of the sum of the weights w(k,n), assuming Gaussian 
independent errors. The error bars therefore indicate the effective weight that each plotted 
point, as the average of several points used in the actual fit, is given to obtain this fit. 

The curves are the result of the model fit, when the model is computed on the same 
values of (k, fi) of all the bins used for the fit, and then averaged in the same way as the 
simulation points for display purposes. This averaging is the reason for the discontinuities in 
these curves in the left panel, which are particularly apparent at low k. The four curves for 
the four values of /x start from the left side at different values of k, depending on the smallest 
k value for which there exists a mode having /x within each bin. The smallest wavenumber, 
with k = k\, exists only for /x = 0 and /x = 1. The model predictions are shown for the 
same values of (fc, /x) as the points, and are plotted as a continuous line only to guide the 
eye. Finally, the cyan and brown dash-dot lines (black for the q 2 = 0 case) are the model 
predictions for /x = 0 and fi = 1 as a function of k (this time, not averaging over any bins), 
shown to indicate the difference with the results in the averaged bins of smallest and largest 
P 

When the points indicating simulation results and the curves showing the model in the 
left panel are divided by the expression [6p^(l + /3/x 2 )] 2 in each of the computational bins, 
using the bias factors obtained in each fit, the non-linear term D(k , /x) is obtained, plotted in 
the right panel. The colored points are now valid only for the free q 2 fit; points for the q 2 = 0 
fit are omitted to avoid excessive cluttering. As before, the points and model predictions are 
averages of the bins used for computing the fit over the coarser bins used to make the plot. 
The averaging of the values and errors is done using the same weights as for the ratio Pf / Pl ■ 

Figure 2 shows that a very good fit is obtained by varying only the 6 parameters in 
equation (3.6), plus the two linear bias factors (colored curves). A value of y 2 = 296.7 is 
obtained, for 552 — 8 = 544 degrees of freedom at redshift 2.2. When fixing x /2 = 0 at z = 2.2, 
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Figure 2. Power spectrum for the 120Mpc//i box at z = 2.2, averaged over the indicated 4 bins 
in /Lt, and with bins in log A: as described in the text. Points with error bars in the left panel are 
the results from the simulation, and colored lines are the fit to equation (3.6) with 6 free parameters 
(in addition to the two bias factors). Black curves are the fit to 5 free non-linear parameters, when 
setting 52 = 0 in D \. The left panel shows the ratio of the transmission power spectrum to the linear 
one, and the right panel shows the non-linear term D(fc,^). Points in the right panel are shown only 
for the fit with 52 as free parameter, to avoid cluttering. Brown dashed curves are the fitted model 
computed at fx = 0 and n = 1, shown to indicate the difference with the averaged bins 0 < fj, < 0.25 
and 0.75 < n < 1; black dashed curves (shown only in the left panel) are the same for the fit with 
<72 = 0 . 




Figure 3. Same as figure 2, but comparing the fit using the D 0 formula (colored curves and points) 
with the one using D\ with free <72 (black curves; black points omitted in right panel), and at z = 3. 


the value of x 2 increases to 349.8. This value of x 2 does not reflect a real “goodness of fit” 
because we have added the parameter e in equation (3.2) to reduce the weight of the modes 
at high k, and because the initial conditions were generated without including the Rayleigh 
distribution in the amplitude of each Fourier mode, so the scatter of the values of Pf/Pl 
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arises only from non-linear coupling. The value of % 2 is therefore substantially less than the 
number of degrees of freedom. Nevertheless, its variation with the model fit can still indicate 
if the improvement of a fit is significant. The fit with q 2 = 0 is already quite good, but it 
improves substantially by including q 2 as free parameter. In particular, the low-A: points are 
better matched for free q 2 , increasing the model power at small k by ~ 15%, implying an 
increase of ~ 7% in the density bias factor. The difference between the two fits for D(k,qi) 
in the right panel goes up to nearly 20% and extends to higher k, because of degeneracy of 
the non-linear parameters with the linear bias factors. 

Results at the higher redshift z = 3 are shown in Figure 3, this time comparing the 
fitting formula D\ with all 6 free parameters (black curves), with the result for Do with 
8 non-linear parameters (colored curves). As before, black points are omitted in the right 
panel. The fits are fairly close to each other for Pp/Pp , with values of x 2 = 252.0 for D i, and 
X 2 = 262.3 for Do, again with a difference at low k where D\ better matches the simulation 
points. In this case, the fit with D\ fixing q 2 = 0 is much closer to being optimal than at 
z = 2.2, with y 2 = 254.6. The fit with the D\ formula is clearly better than with Do, with 
two fewer parameters; we have found this to be true also at z = 2.2. Moreover, we see in the 
right panel that the formula Do converges very slowly toward unity at small k. The reason 
is that the value of the a n i parameter in equation (3.4) obtained for this fit at z = 3 is small, 
a n i ~ 0.2, and there is a strong degeneracy with the value of the linear bias factors. This 
formula therefore easily leads to unphysical parameter values when the linear bias factors are 
not determined independently from the fit to the power spectrum of a simulation of limited 
box size. The values of the parameters are listed in the tables in Appendix B. 

4.1 Results for the linear bias factors 

In the limit of small k, the ratio Pf/Pl approaches b 2 F Al+/3qi 2 ) 2 . This is indeed the behavior 
shown by our results in the left panel of Figures 2 and 3. The size of the box limits the number 
of Fourier modes available in the simulation at low k and hence the accuracy to which the 
linear bias factors can be measured from our fit. In addition, degeneracies between these 
bias factors and the non-linear parameters are present when fitting the simulation results for 
Pf/Pl- 

The values of the bias factors and the redshift distortion parameter are shown in Figure 
4 at the 5 redshift outputs that we will use in most of our models, and for the three cases 
we have considered: the D\ formula with free q 2 , fixing q 2 = 0 in D \, and the Dq formula. 
The two physical bias factors, b T $ and b Trj , are obtained from the transmission bias factors 
bps and bp v derived from our fits of Pp(k, qi)/Pp(k) through equation (2.4), and are shown 
in the center and right panels. Errorbars are derived from the Monte Carlo Markov Chain 
computed for the fits, and they generally overestimate the purely statistical errors for a fixed 
fitting function because our y 2 function is too low due to the e parameter and the absence 
of Rayleigh-distributed amplitudes in the simulation. However, they are dependent on the 
fit model and they do not include the systematic errors that are investigated in §5. 

The results of the bias factors for Do differ substantially from the ones derived with 
D i, and they are less reliable because the function D{k,[P) converges too slowly towards 
unity at small k, as seen in Figure 3. This slow convergence is unphysical because non-linear 
effects are actually very small on the largest modes of the L120 simulation box. This problem 
for Do is worse at high redshift and, together with strong degeneracies with the non-linear 
parameters, is the cause of the strange redshift dependence of the bias factors. For this 
reason, we focus in the rest of this paper on results obtained with D\. 
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Figure 4. Results of the fitted values for the two bias factors b T $ and b TTJ , and the redshift distortion 
parameter /3 as a function of redshift, for the 5 redshifts outputs of the L120 simulation. The error 
bars are ler, as returned from the MCMC fitting. To avoid superposition of errorbars, green points are 
shifted slightly left and blue points slightly right, but they are all at our 5 standard redshift outputs. 
Results are shown for the fitting formula D\ with all 6 free non-linear parameters (red squares), for 
D\ setting 92 = 0 (blue circles), and for the formula Dq with 8 free parameters (green triangles). 
The differences among the fitting models are caused by degeneracies with the non-linear parameters, 
which are more severe for Dq. In the right panel, the absolute value of the radiation bias parameter 
b T r is also shown as crosses. 


The redshift distortion parameter (3, shown in the left panel, is predicted to have a 
value near 1.4 at the most commonly observed redshift in BOSS, z — 2.3, and to decline 
with redshift. There is little variation of the value of [3 when fixing <72 to zero. The value we 
predict is slightly smaller than that of M03, who found /3 = 1.58 at z = 2.25. Our result for 
(3 for the Dq formula is closer to that of M03, but this is not for the same reason since an 
independent method was used by M03 to determine the bias factors. The dependence of (3 
on the fitting method and physical parameters will be further discussed in §6. 

The physical bias parameters have a relatively weak dependence on redshift. Previous 
results reported in terms of the transmission bias factors indicated a very rapid evolution 
with redshift, due mostly to the change in the mean transmission F. The predicted value 
of b T s — 0.6, nearly constant with redshift, has the physical meaning that the Lyct effective 
optical depth fluctuates on large scales by ~ 60% of 5, the fluctuation in the mass density 
field. As mentioned above, the density bias factor drops by ~ 8 % at z = 2.2 when fixing 
52 = 0, corresponding to the ~ 15% variation of power at low k seen in Figure 2. This is an 
indication of the uncertainty due to the degeneracy with non-linear parameters and the use 
of different fitting formulae. 

The bias factor of the peculiar velocity gradient is computed from equation (2.6). For 
the logarithmic derivative of the growth factor, we use the values for the cosmological model 
of our fiducial simulation: f(tl) = (0.9875, 0.9895, 0.9911, 0.9924, 0.9935) at z = (2.2, 2.4, 
2.6, 2.8, 3). The value of b rn is below unity and decreases with redshift. This means that the 
Lya forest behaves differently from a large-scale structure survey of objects with a selection 
function that is independent of the peculiar velocity gradient along the LOS, r /: the effective 
optical depth fluctuates only by ~ 70% of the fluctuation in 77 at z = 3. We shall return in 
§7 to the physical reason why b rrj is less than unity and decreases with redshift, which is a 
general characteristic in the results of all our simulations. 
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Finally, we compute also the radiation bias factor b T r, defined at the end of §2, using the 
probability distribution of F in the simulation LI20 at each redshift output, and equation 
(2.8). The results are shown as crosses in the right panel of figure 4. The values are a 
factor of ~ 2 smaller than b rr] . This shows that the approximation proposed by [51], in 
which b Tri = b T r, fails in an important way. Following on our discussion at the end of §2, 
we believe the reason is that non-linear evolution of the small scale fluctuations in the Lya 
forest substantially modifies the value of b Tr) . 

As far as the linear power spectrum is concerned, the main goal of numerical simulations 
of the Lya forest should be to accurately predict the value of (3(z) and b T1J (z ), and to examine 
the model dependence of these functions, to compare to observational determinations. Our 
results in this section already show the main difficulty involved in this goal: the numerical fits 
we obtain depend on the fitting formula that is used, and on the value of the e parameter we 
have chosen in equation (3.2) to increase the importance of the low-/c modes. When setting 
<72 = 0, the fit to the low-A: points is worse and both b T s and b TT1 go down. However, if q -2 
is set free then the q\ parameter determining the limit of the non-linear correction at low 
k becomes highly degenerate. Determining accurate values of the bias parameters with a 
reliable non-linear correction can only be done with large numbers of simulations on large 
boxes, to have better statistics for the power at low -k. Since we have only one simulation 
with L = 120 h Mpc and a few on smaller boxes, the results in this paper on linear bias 
factors are of limited accuracy, as exemplified by the difference in our two different fits with 
D x . 



Figure 5. Results of the non-linear parameters of D\ as a function of redshift, for the 5 redshifts 
outputs of the L120 simulation. The error bars represent the ler contour as obtained from the MCMC 
fitting. Solid red lines are for all 6 parameters left free, and dashed blue ones are for fixing q 2 = 0. 
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4.2 The non-linear part of the power spectrum 

The values of the non-linear parameters are shown for the fitting model D\ in Figure 5, for 
all 6 parameters being free (red solid line), and for = 0 (blue dashed line). We plot 
instead of k v to reduce the amount of degeneracy among parameters. We comment here on 
the principal features of the results on these parameters and what this reflects on the shape 
of the non-linear function D(k, fi). 

The scale k p is generally near 10/i/Mpc, or k p /H ~ 0.1 s/km. This is a characteristic 
value for the Jeans scale of the photoionized gas in the IGM, as discussed in §3.4. This 
parameter is used to fit the declining power with k that occurs above this scale even for 
[i = 0, clearly seen in Figures 2 and 3. The value of q±, controlling the amplitude of the 
non-linear power enhancement, is q± — 0.6 when no second-order term is included, but a 
substantial degeneracy occurs with the parameters q 2 , k v , k p and a v when q 2 is included in 
the fit, causing large changes of q\ that vary with redshift. 

Non-linear effects imply a change of the sign of the quadrupole of Pp(&, /^) as k increases. 
At small k, the power is largest at fi = 1 due to the Kaiser effect, and at large k the effect of 
velocity dispersion takes over and the power is largest at fi = 0. The quadrupole is zero at a 
point where curves of different p, cross each other. This point reflects a characteristic scale at 
which non-linearity changes the sign of the power spectrum anisotropy, and shifts to the left 
(larger scales) as the amplitude of the mass power spectrum increases from z = 3 to z = 2.2. 




Figure 6. Left panel: Pp(/c,/x)/P/,(fc) as a function of p, 2 , for the bins centered at k = 
(0.56,1.44, 3.08, 6.58) /i/Mpc (or bins number 1, 6, 10 and 14 above k t ). Right panel: same points for 
D(k,n) as a function of /x 1 ' 5 . 

The dependence of the non-linear correction on fj, is very well represented by the power- 
law //'" inside the exponential in equations (3.4) and (3.6), where b v is close to 1.5 at all 
redshifts and is subject to little degeneracy with other parameters. To see in greater detail 
the dependence of the power on /x, Figure 6 shows Pf/Pl as a function of ^x 2 in the left panel, 
and D(k,fi) as a function of /x L5 in the right panel, for a set of four selected bins in log A:, 
corresponding to bins number 3, 7, 11 and 16 in Figure 2 starting from the right edge. The 
four values of k of these bins are indicated in the figure, in units of h/ Mpc. This figure shows 
the actual value of Pf that we use in our fits for the 16 bins in fj, in our analysis, instead of 
the averages over coarse bins as in previous figures. The results are shown at z = 2.2 and for 
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the D\ fit including q 2 as free parameter. We see in the left panel that at the scale where 
the quadrupole changes sign (k = 3.04/i/Mpc), the power is indeed almost independent of 
p. The right panel shows that the ^-dependence of the non-linear correction to the power is 
surprisingly well fit by a simple power-law at all values of k. 

It will be interesting to express our non-linear parameters in terms of the characteristic 
scale where the anisotropy changes sign. We define the wavenumber k na (the subindex is for 
the scale of non-linear anisotropy ) as the one that obeys 


Ppiknai H — 0 ) — Ppiknaik 1 — 1 ) ■ 

The /a dependence of the transmission power spectrum for our fitting formula D\ is 

Pp(k,n ) oc (1 + Pp 2 ) 2 exp -p a (k)p b 

where, if we use for simplicity the formula D\ with q 2 = 0, we have 

/ k \ av 

Pa(k) = qiA 2 (k) 

The value of k na can be computed from our non-linear parameters from equation 

Pa {kna) = 21og(l + /3) . 


(4.1) 


(4.2) 


(4.3) 


(4.4) 


As it turns out, the function of p in equation (4.2) happens to be nearly constant in the range 
0 < /x < 1 when b v ~ 1.55 and /3 ~ 1.4, if p a obeys equation (4.4). As we shall see, these 
values of b v and P do not change much for all the other models we examine in this paper. 
This is the reason that all the curves for Pp/Pl with different p cross each other nearly at 
the same wavenumber k na for all our models. 

We compute the value of k na for all our models with q 2 = 0, by iteratively solving the 
equation 


k 


na 


(27r) 2 k“ v log(l + P) 
qi PLpkna) 


1/(3+a-u) 


(4.5) 


The parameter values of our fits are given in the tables of Appendix B. The parameter k na 
is particularly useful because its correspondence with the crossing point of the curves means 
that it has very little degeneracy with all other parameters, and therefore its error is small. 

We have not been able to think of an analytical explanation for the surprisingly well 
matched /i-dependence of log D(k,p) to the power-law p bv , with b v ~ 1.55. We note that 
b v <2 implies a singular second derivative of Pp(k,p) at p, = 0. 

Finally, we also note that once q 2 is introduced as a free parameter, the power-law index 
a v becomes very small, particularly at low redshift, and a large degeneracy is introduced 
among the parameters k v , a v , q\ and q 2 ■ The improvement in the fit obtained by including q 2 
is modest, and this improvement is further reduced for simulations in smaller boxes, where 
less information is available on the /c-dependence (and the parameter degeneracy increases). 
We therefore will fix q 2 = 0 for our fits to the results of most of our simulations. 


5 Convergence Tests: Resolution, Box Size and Numerical Method 

This section addresses the degree to which our fit results from the various simulations we 
analyze have converged when box size, resolution, grid size and numerical method are varied. 
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Variations of the power spectrum with physical parameters of the Lya forest will be discussed 
in §6. 

We remark here one important difference between our work and that of M03. We do 
not use any splicing technique (see also [7]) to combine results from simulations of different 
resolution or box size, a method introduced by M03 to attempt to better reach a convergent 
solution. In this work we simply fit both the linear bias factors and the non-linear parameters 
of the analytic formulae of equations (3.4) and (3.6) to the power spectrum of each simulation. 
Our results can probably be improved by using this splicing technique, in particular for the 
linear bias factors, but we did not have enough simulations in this work to do this for all the 
models we analyze, and we left this for future studies. The emphasis of our analysis is more 
focused on the variations of the bias parameters and the non-linear form of the Lya power 
spectrum with respect to a reference model (our fiducial simulation), rather than aiming for 
a highly accurate convergence of the results, for which a larger number of simulations on 
large boxes are required. 

Before we start, we also note that all our SPH simulations generally use the same 
random initial conditions for the phases of Fourier modes of fixed kL. Moreover, as described 
in §3.1.1, the amplitudes of the modes are always equal to the variance predicted by the 
power spectrum, instead of being generated randomly with the Rayleigh distribution. This 
minimizes the random sampling variance due to the finite box size and allows for a more 
direct comparison of different models. In general, however, there is no guarantee that the 
average Pp obtained by suppressing the Rayleigh distribution of amplitudes in the initial 
conditions is the same as the correct Pp that can be derived only by including the Rayleigh 
distribution, and this will need to be tested in future studies. As we shall see in this section, 
the variations of the low-A; power spectrum from simulations of the same model but different 
random initial conditions are similar to the differences introduced by resolution and box size, 
and the overall accuracy to which we can measure the bias factors from our simulations is 
not better than ~ 10%. 

5.1 Resolution 

We start by checking the dependence of our results on the resolution. There are three different 
quantities in relation to the resolution of the simulation and the analysis that is done for 
computing the power spectrum that must be tested for convergence: the number of particles 
in the simulation, the number of cells in the grid to compute the hydrodynamic variables, 
and the number of pixels in the simulated Lya spectra. 

We compare first results for the power spectrum from the same simulation and spatial 
grid, but varying the pixel resolution of the computed Lya spectra. Results for Pp/Pp are 
shown in figure 7, where black curves are for our fiducial model and colored curves and points 
are for the P1024 model. The fiducial model is the same as the one used in §4, but with the 
box size L = 60 h~ 1 Mpc, and the P1024 model is exactly the same but with the number of 
pixels in the Lya spectra doubled to 1024 (Table 3.1). Black points are omitted to avoid 
cluttering, but they follow the black curves similarly to the colored ones. The P1024 model is 
computed by calculating the optical depth in 1024 pixels along the LOS, and then averaging 
exp(—r) for every two pixels to obtain a cubic grid of 512 3 values of the transmission fraction, 
on which a Fast Fourier Transform is done in the same way as for all our other models. Both 
models are fitted with our D\ formula with q 2 = 0. 

The results are shown at z = 2.6, and are very similar at other redshifts. The increase in 
the pixel resolution practically does not affect the transmission at the level that is discerned 
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Figure 7. Comparison of the power spectrum result obtained for our fiducial simulation and model 
analysis (black curves; box size L = 60 h _1 Mpc, 512 3 particles and grid cells, and 512 spectral pixels), 
with the case of doubling the spectral pixels to 1024, at redshift z = 2.6 (colored curves). Black points 
are omitted to avoid cluttering. 


in figure 7 (the variation is less than 2%). Clearly, the spectral pixel resolution is sufficient 
for our purpose, and any errors it is causing are much smaller than those due to the grid size 
and number of SPH particles in the simulation. 

Next, we compare two analyses of a simulation using spatial grids of different resolution. 
For that purpose, we take our highest physical resolution simulation, R640 (equal to the 
fiducial one but with 640 3 particles instead of 512 3 , although still analyzed on a 512 3 cells 
grid; see Table 3.1), and compare it to R640C, the same simulation analyzed on a grid of 
only half the sampling, with 256 3 cells. The results for Pf/Pl are shown in figure 8, at 
redshift z = 2.2. Larger differences can be appreciated in this figure. First of all, a coarse 
grid increases the power on small scales, k > 5/i/Mpc, probably due to the noise introduced 
by a grid cell of 0.23 h _1 Mpc in R640C. This also generates a roughly constant decrease of 
the power on large scales in R640C compared to R640. In fact, for k < 0.5 h~ l Mpc, the black 
points are systematically higher than the colored ones by ~ 7% at all values of fi. The change 
in the small-scale power alters non-linear couplings to produce an impact on the large-scale 
power that becomes constant on very large scales, therefore altering the bias factors. Results 
at other redshifts are qualitatively similar. 

Figure 8 also exemplifies the uncertainty in our fits for recovering the low-A: power: the 
colored curves move above the black ones at high fi and low k, which is opposite to the 
simulation results for Pf/Pl indicated by the points. This is again reflecting the uncertainty 
due to the fitting formula that is chosen. The fitted curves can better reflect the behavior 
of the simulation points when the q 2 parameter is left free, but the reliability of the derived 
bias factors is still subject to a similar error because the small box size does not allow one 
to obtain the bias factors more accurately. 
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Figure 8. Comparison of Pf/Pl at z = 2.2, for the R640 simulation with 640 3 particles analyzed 
with a grid of 512 3 cells (black curves and points), with the same simulation analyzed with a coarser 
grid of 256 3 cells (R640C, colored curves and points). Errorbars for the black points are suppressed 
to avoid cluttering, but they are identical to the colored ones. 




Figure 9. Redshift distortion parameter and bias factors fitted for three simulations of increasing 
physical resolution: R384, fiducial, and R640. These three models all have a 512 3 grid, and we 
compare them also to R640C, with a coarser grid of 256 3 and the same particle number as R640. 
Blue circles for R384 are shifted slightly left, and green triangles for R640 slightly right, to avoid 
superposition of errorbars. 


The bias factors derived from the q 2 = 0 fits are shown in figure 9. The errors are the 
formal ones derived from the x 2 fits, but the more important errors on the biases are of a 
systematic nature arising from the limited size of our boxes, and the need to fit the curves of 
Pf/Pl at small scales where the simulations provide much more information than at large 
scales. The results are fairly similar for all the models shown in this figure. The largest 
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difference is the higher value of f3 derived for R640C, but as discussed above, this does not 
reflect the true difference of power at low k between R640C and R640 shown by the colored 
and black points in figure 8, and is an artifact of our fit with q 2 = 0. The real systematic 
errors in the bias factors are better reflected by the difference with the fits with free q 2 - 

We conclude that the errors introduced by the limited resolution of the grid size are 
substantial, and that future work should study the convergence with both particles and grid 
size for safely obtaining results for Pp to accuracies better than 10%. We note that our L80 
and L120 simulations are analyzed also with a 512 3 grid, which for a larger box size implies 
a lower physical resolution, with an expected impact on the power spectrum and bias factors 
following that seen in figures 8 and 9. 




Figure 10. Comparing the power spectrum of simulations R640 (colored curves and points) and 
R384 (black curves, points omitted), at z = 2.6, with higher and lower resolution compared to the 
fiducial model, respectively. Left panel: fits with 52 — 0. Right panel: fits with free < 72 - 

Finally, we compare simulations of different resolution, varying the number of SPH 
particles. The simulations R640 and R384, which are the same as the fiducial one but 
varying the number of particles to 640 3 and 384 s , respectively, produce results for Pp/Pp 
that are shown in figure 10 at z = 2.6. In general, the changes due to the physical resolution 
are comparable to those caused by the resolution of the spatial grid but of opposite sign: 
decreasing the particle resolution decreases the power at high k, although in this case this 
is clear only for low /i. In addition, the power at low k generally drops with increasing 
resolution. In this case, we show in the left panel the two fits with q 2 = 0, and in the right 
panel the two fits with free < 72 - The fit at low k follows the points better for free q 2 , as 
expected because of the extra degree of freedom, and the difference between the curves also 
better reflects the differences in the points for R640 and R384. However, as discussed before, 
it is not clear from the present simulations if the fits with free (72 (which have a lower value 
of q\ and therefore converge faster to the linear solution) are a more reliable result for the 
bias factors and the non-linear correction D(k,p) at low k. 

The derived bias factors for q 2 = 0 are also shown for both models in figure 9. The bias 
factors for free q 2 are not shown, but are similar to those of the fiducial model shown below in 
figure 14. In general the density bias factor decreases by ~ 5% when the particle resolution 
is doubled at fixed grid size, and increases by a smaller amount when the grid resolution is 
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doubled. We note here that the changes at low k induced by the small-scale resolution effects 
are due to non-linear couplings that depend on the fact that we always keep a fixed value 
of F ; the result would in general be different if instead we kept a fixed value of the ionizing 
background intensity. 

We therefore conclude that the limited resolution and grid size in our simulations affect 
the values of the bias parameters by ~ 5%, and unfortunately our results do not show yet 
a clear convergence with the resolution. Increasing the number of particles or the number 
of grid cells modifies the results in opposite directions. Higher resolution simulations will 
be needed to better understand the conditions for convergence and to obtain more accurate 
results for the non-linear power Pp(fc, //) and the bias factors. Nevertheless, our derived values 
of the bias factors and ft at any redshifts vary by less than 10 % in all of our simulations, and 
the predicted non-linear shape is also not subject to greater relative changes arising from 
resolution, so we believe that our results are reliable within this level of uncertainty. 

Our results on the sensitivity to the simulation resolution can be compared to those 
obtained in [35]: their figure 12 shows the quantity A 2 F = fc 3 Pp/(27 t 2 ) (equal to what we plot 
below in figure 19). In agreement with our results in figure 10, lower resolution decreases 
the power at high k and increases it at low k, and these changes are more pronounced at 
low /z. However, the sensitivity to resolution seems to be weaker in our SPH simulations. 
Our R640 and R384 simulations have an interparticle separation that is about twice the cell 
size of the T10 at 128 and L10 at 256 simulations of [35], and yet the differences between the 
two simulations are ~ 40% in [35], and only ~ 10% between R640 and R384 in our case. 
The likely explanation for this is that the value of ~ F is not held fixed when comparing 
simulations of different resolution in [35], and that our results shown at z = 2.6 are mostly 
sensitive to gas densities above the mean density of the universe, where the resolution of 
particle-based codes improves over that of Eulerian codes. At higher redshift, however, the 
Lyct forest probes lower overdensities and the resolution effects in our simulations probably 
become more severe. 

5.2 Simulation Variance and Box Size 

The limited box size of any cosmological simulation implies the presence of a statistical error 
and a systematic error on any quantities that are measured from the simulation, such as 
the power spectrum of any tracer. The statistical error is due to the shot noise arising 
from the finite number of structures of a given scale that are formed within the simulation 
volume; if a number N s i m of simulations with the same box size are performed with different 
random initial conditions, this statistical error is reduced as N ■ ' . The systematic error, 
which is not reduced by averaging over several simulations, is due to the discretization of 
Fourier modes, which need to have cartesian components on the box axes that are multiples 
of k\ = 27 t/L, and the absence of any modes below the wavenumber k\. We therefore start 
by looking at the statistical variations of the power spectrum in the fiducial model, with 
L = 60 1 Mpc, by examining three simulations run with three different random seeds for 

the initial conditions: Seed 1 is for our fiducial model shown in previous figures, Seed 2 is 
for independent initial conditions where the Fourier mode amplitudes are generated with the 
correct Rayleigh distribution, and Seed 3 for independent initial conditions varying only the 
Fourier mode phases of Seed 1, but keeping the square amplitudes fixed to the mean value 
determined by the matter power spectrum, as in Seed 1 (see §3.1.1). 

Figure 11 shows Pf/Pl at z = 3 for the first two seeds as points, and the fits with (72 = 0 
as curves. Seed 2, with the Rayleigh distribution, is shown as colored curves and points, and 
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Figure 11. Comparison of Pf/Pl for two different random initial conditions for the fiducial simula¬ 
tion at z=3, with a 60 Mpc /h box size. Seed 1 is our usual fiducial simulation (black curves and points 
with errorbars omitted), with all Fourier mode square amplitudes fixed to the mean value predicted 
by the matter power spectrum of the model. Seed 2 is for different initial conditions with the Rayleigh 
distribution of amplitudes (colored curves and points with errorbars). 





Figure 12. Comparison of the redshift distortion parameter and linear bias factors for three different 
random initial conditions for the fiducial simulation, obtained from the fits with 52 = 0. 


our usual fiducial simulation (Seed 1) is the black curves and points. The difference in the 
two fits is very small, the largest difference occurs at low k and high ^ and is only ~ 3%. 
This is despite the large random differences seen at low k due to the introduction of random 
amplitudes for seed 2. The fit does not have much freedom to vary owing to the large number 
of modes at small scales, resulting in very small variations for the high-A; points between the 
two simulations, so the large scatter at low k gives small variations in the derived bias factors. 
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Despite this, we see that the fits are generally good for both seeds. At 2 = 3, the (72 = 0 
fits are in fact better at low-A; than the cases at lower redshift shown in previous figures, and 
are also closer to the free <72 fits. Note that the even though the errorbars are shown around 
the measured simulation points, they need to be computed for the predicted Pp of the fitted 
model, so points that are below the model are not as statistically discrepant as they look in 
the figure (e.g., fourth red circle and first green triangle starting from the left). 

The derived linear bias factors are plotted in Figure 12 at all redshifts, for all three 
seeds, for the (72 = 0 fits. Typically, the two bias factors can have random variations of up 
to ~ 4% when varying the initial conditions, and this is not strongly affected by including 
the Rayleigh distribution of Fourier mode amplitudes. We have tested that these conclusions 
remain valid for an additional two independent seeds, not shown here. This uncertainty due 
to the random initial conditions is therefore comparable to the one due to the resolution 
of the simulations. These tests also show that the difference between including or not the 
Rayleigh distribution of mode amplitudes in the derived average Pp is small enough for the 
purpose of the present study. 

Next, we compare two simulations of the same model that differ only in the box size, 
but have the same physical resolution in both particles and grid cells: L120 and R384C. The 
R384C simulation has box size L = 60 /r _1 Mpc, half of that of L120, and also half the number 
of particles and grid cells along a LOS, so that the initial interparticle spacing and grid cell 
size are the same. The ratio PpjPp is shown in figure 13 for these two models at z = 2.2. 
The upper panel shows the fits with free ( 72 , and the lower panel the fits with q 2 = 0. Two 
types of differences induced by the finite box size are seen in this plot. First, at low k, the 
accuracy in fitting bias factors is increased as the box size is increased and more points are 
available to determine the linear limit of the power spectrum. This is seen more clearly for 
the lowest ^1 curve, which is raised by ~ 10% from R384C to L120 in the free q 2 fit. Second, 
the points for the R384C simulation are generally lower than for L120 even at high k, with 
the largest difference ocurring also for the low fj, curve, where it reaches ~ 5%. The main 
reason for this difference at high k is that we are always computing the Lya transmission 
power after fixing F. and the value by which we need to multiply the optical depth to achieve 
a fixed F changes in the various models we show. The large-scale power that is missing in 
the R384C simulation results in larger voids with lower densities, and for a fixed intensity 
of the ionizing background, this would result in a higher F. The intensity of the ionizing 
background therefore needs to be lowered to have the same F in the two simulations, and 
this induces the change in power at high k. At the highest values of k near 10 h —1 Mpc, the 
power measured in the two simulations shown by the colored and black points is very close, 
but as we have seen previously this power is quite sensitive to the particle and grid resolution, 
which are the same in these two simulations. 

Linear bias factors are shown in figure 14 for the fits with free q 2 (blue circles for L120 
and purple rhombi for R384C), where a third model of intermediate box size, L80, has been 
added (green triangles). The results for the fiducial model, on a 60/i _1 Mpc box, are also 
shown as red squares (these were shown only for the (72 = 0 fit in previous figures). In 
general, the bias factors are higher by ~ 10 % for the free (72 fit compared to the (72 = 0 fit. 
While the bias factors are nearly the same for LI20 and R384C, the L80 simulation has a 
larger value of b T s by ~ 5% at all redshifts, while b TTI is nearly the same, and therefore (3 is 
lower for L80. This may seem strange because L80 is a simulation of intermediate box size. 
The L80 simulation also has the same particle resolution as the others, but its grid size is 
smaller (with the same number 512 3 of cells as the larger box of L120). As seen previously 
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Figure 13. Comparison of Pp/Pp from simulations L120, with box size of 120 Mpc/h (colored points 
and curves), and R384C, on a 60 Mpc /h box (black curves and points with omitted errorbars), at 
z = 2.2. The two simulations have equivalent physical resolution, with 768 3 and 384 3 particles, and 
512 3 and 256 3 grid cells, respectively, at 2 = 2.2, so differences should be due to the effects of box 
size and the random initial conditions displaced by a factor of 2 in k only. Upper panel shows the fit 
with q 2 as a free parameter, and bottom panel is for <72 = 0. 


in figure 9, increased grid resolution raises the value of b T s and lowers (3, and correcting for 
this brings the green lines (for L80) a bit closer to the blue and purple lines (for LI20 and 


- 30 - 


























































Figure 14. Comparison of the bias and /3 parameters for different box sizes. Notice that while the 
R384C and L120 simulations have equivalent physical resolution, L80 has the same particle density 
but a different physical grid size. The results shown here are for free <72 in the D\ formula. 


R384C) in figure 14. However, most of the difference remains, and we believe this can be 
assigned to the effects of random initial conditions. Even though these simulations have the 
same Rayleigh-suppressed initial conditions on modes of fixed kL, the power of a simulation 
at fixed k still has random variations, and therefore we can expect some of the variations in 
bias factors analogous to the ones between seeds 1 and 3 seen in figure 12. 

We conclude that the systematic effects of the box size are comparable to the random 
variations due to different initial conditions seen in figure 12, and also to the systematics due 
to limited resolution in our simulations. If anything, the error resulting from using different 
fits to the whole shape of Pf/Pl , which we have illustrated by including or not including the 
additional free parameter t /2 in this paper, is a bit larger: both b T s and b Trj are systematically 
larger by ~ 10% for free (72 than for fixed <72 = 0 at low redshift in most of our models, with 
a smaller difference at high redshift. The better fit that is obtained for free c /2 to the low -k 
points suggests this to be more reliable, but this fit also implies a rather small value of the 
qi parameter (see figure 5), and therefore a surprisingly fast convergence of the non-linear 
factor D\ to unity at low k. This remains to be tested with more large box simulations. 

5.3 Comparing Lagrangian and Eulerian simulations 

As a final test, we analyze here a simulation that is directly run on an Eulerian grid instead 
of using gas particles. This is the simulation labeled ” Euler” in table 3.1 and described in 
section 3.1.2. As a comparison, we have run an SPH simulation designed to match the physical 
parameters of the Euler simulation, so that a comparison can be made. This simulation is 
labeled ”Lagrange” in table 3.1: it has the same box size of L = 50 h^ 1 Mpc, and resolution of 
512 3 both in number of particles and grid cells. The model mass fluctuation power spectrum 
is the same, and we have also attempted to match the two simulations to have the same 
density-temperature relation, which depends on the model of the ionizing background and 
the heating and reionization history. Unfortunately, it is not possible to do this match exactly 
owing to the different methods for computing heating in the two simulations. In addition, 
the initial conditions were not the same in the two simulations, (in this case both of them 
include the Rayleigh distribution of Fourier mode amplitudes), so this inevitably introduces 
some difference in the results. 
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Figure 15. Mean temperature as a function of density for the Eulerian and Lagrangian simulations, 
at 2 = 2.3 and 2 = 2.9. 

The values of the mean transmission for the Eulerian simulation were determined by 
the model of the ionizing background evolution used in that simulation, and are equal to the 
following values: F = 0.8042, 0.7484 and 0.6725 at the three redshift outputs of z = 2.3, 2.6 
and 2.9, respectively. We choose to rescale the optical depths in the Lagrange simulation to 
match these same values of the mean transmission, which are different from our standard ones 
used in previous figures. The redshift outputs for the simulation Lagrange are our standard 
set of 5 values, for which we set the mean transmission to the same value as the Eulerian one 
at z = 2.6, and to values following a power-law dependence of the same form as in equation 
3.1, separately chosen to match the mean transmission of the Euler simulation for z < 2.6 and 
for z > 2.6. This results in the following values: F = (0.8212,0.7864,0.7484,0.6989,0.6452) 
at z = (2.2, 2.4, 2.6, 2.8, 3.0). 

The full density-temperature relation of the Euler simulation is shown in figure 16, at 
redshifts z = 2.3 and z = 2.9. The temperature that is shown is the average in all the cells 
that have the density in each bin of width A log p = 0.1. The relation is not exactly a power- 
law, and it steepens above p/p ~ 1 due to shock-heating of the gas in the sheets and filaments 
of the Lya forest as structure formation develops. The two simulations are fairly close to 
each other, but in the low-density regime the Eulerian one has lower temperature by ~ 0.1 
dex. At densities p/p > 1 the shock heating in the Eulerian simulation is higher than in the 
Lagrangian one, particularly at low redshifts, and at p/p > 5 the Lagrange simulation shows 
more heating at low redshift. This seems to be affected by the collapse and sudden heating 
of a few clusters in the simulation and is therefore affected by random variance depending 
on initial conditions. We have obtained a linear regression fit to the log T — log p relation 
of the Euler simulation in the low-density regime —1.5 < log (p/p) < —0.3, in which shock 
heating is not important, with the result indicated in Table 3.1 for z = 2.6, which has a 
weak variation with redshift. The difficulty in obtaining a more precise match of the density- 
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Figure 16. Comparison of Pf/Pl for the Euler and Lagrange simulations, using an Eulerian code 
and the GADGET-II code, of the same cosmological model, at z = 2.6. 





Figure 17. Redshift evolution of (3 and the bias factors for the Euler and Lagrange simulations, from 
the fits with 92 = 0 . 


temperature relation in the Eulerian and Lagrangian simulations illustrates the problems that 
are encountered when testing that numerical simulations using different codes and different 
methods make the same predictions for the Lya forest power spectrum. More precise tests, 
that start from exactly the same initial conditions and force the same thermal histories, 
should be performed in order to investigate the differences between the two frameworks more 
completely (see, e.g., [50], where a comparison between SPH and Eulerian codes in terms of 
the ID transmission power spectrum is presented and where differences of order < 5% are 
found). 

The ratio Pf/Pl is shown in figure 16, for the Euler simulation as colored curves and 
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points, and the Lagrange one for black curves and points with omitted errorbars, at z = 2.6. 
The curves are fits to the D\ formula with <72 = 0. Important differences between the two 
simulations are apparent: first, at low k, the high-// curves are higher for the Euler simulation, 
implying a higher predicted value of /3. At high k, the curve at lowest n is lower for Euler 
than Lagrange, whereas the high-// curves are also at lower power except at k > 5 h _1 Mpc, 
where the power decreases more slowly for Eulerian simulation. The differences are clear and 
at the level of ~ 10 %, and they may be due to several effects having to do with the different 
resolution and methods of the simulations. It is useful in particular to compare with figure 
10 , which shows that poor resolution produces a sharper decline of the power at the highest 
k, suggesting that the poorer resolution of the Lagrange simulation compared to the Euler 
one at mean density can explain this difference in figure 16. Nevertheless, the fact that the 
two simulations agree within 10 % on the predicted power (except at low k and high /r where 
the difference grows to 20 %, but is sensitive to random variance from the different initial 
conditions) is reassuring. 

The bias factors are shown in figure 17. The redshift distortion factor (left panel) pre¬ 
dicted by the Euler simulation is substantially higher than for the other simulations analyzed 
so far, corresponding to the higher power of the colored curve at low k and high /r in figure 
16, whereas the variations in b T s are very small, but as discussed above this is sensitive to the 
variance in the different initial conditions. The differences in the predictions between the two 
methods will need to be understood in more detail before theoretical predictions can be made 
from these hydrodynamic cosmological simulations at a high level of accuracy. However, the 
main features of the evolution of these bias factors are common to all our results: b T s is 
roughly constant with redshift and near 60%, and b rri is slightly below unity and declines 
with redshift, implying also a decline of /3 with redshift. 

6 Dependence of the Lya Power Spectrum on the Physical Model 

We analyze in this section the dependence of the Lya transmission power spectrum on the 
physical characteristics of the IGM and the cosmological model. As described in §1, the 
Lya power spectrum depends mainly on the amplitude and shape of the initial matter power 
spectrum, the density-temperature relation, and the value of the mean transmission at each 
redshift. Other properties of the cosmological model should not affect Pp(/c,/r), except for 
rescalings: for example, changes in Hq, and simply imply a rescaling of the Jeans scale 
and any physical scale in the power spectrum in terms of angular and redshift separations, 
which are the observed coordinates of spectral pixels. The details of the gas cooling time, 
which depend on fand the ionizing background spectrum, can be thought of as being 
incorporated in the density-temperature relation and its history, which determine the impact 
of the Jeans scale and the thermal broadening on the Lya power spectrum. 

The dependence of the Lya power spectrum with the most relevant physical parameters 
are described below, by examining the following cases: 

• Dependence on the amplitude of the power spectrum, characterized by the parameter 
erg, which we vary for our fiducial model (§ 6 . 1 ). 

• Dependence on the power spectrum slope on the range of scales examined by our 
simulations, which we illustrate by comparing the fiducial and Planck models (§6.2). 

• Dependence on F at a fixed redshift, and redshift evolution (§6.3). 
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• Dependence on the density-temperature 
(§6.4). 

Throughout this section, we use only fits 

6.1 Power spectrum amplitude 



relation, parameterised as T = To(pfp) 1 1 
the formula D\ fixing t /2 = 0. 



Figure 18. Comparison of Pf/Pl (left panel) for models with different power spectrum amplitude, 
SO.76 with a 8 = 0.7581 (colored curves and points), and the fiducial model with a s = 0.8778 (black 
curves; points omitted), at z=2.6. The non-linear term D(k,p) is shown in the right panel. As ex¬ 
pected, D — l increases with the amplitude, and the characteristic scale at which non-linear anisotropy 
starts dominating, where curves with different p cross each other in the left panel, increases ( k na de¬ 
creases) with the amplitude. 

The ratio Pf/Pl is compared in figure 18 (left panel) for the model SO.76 and the 
fiducial model. The only difference between these two models is in the amplitude of their 
linear power spectrum, which is lower by a factor (0.7581/0.8778) 2 = 0.746 in SO.76. The 
colored curves are generally higher than the black ones by a factor that is slightly smaller 
than 0.746, indicating that Pf increases slowly with erg, so Pf/Pl decreases although not as 
fast as ag 2 . As the power spectrum amplitude increases, we expect that non-linear effects 
become important at increasingly long scales. This is confirmed by the fact that the char¬ 
acteristic wavenumber k na at which the curves with different p cross each other decreases 
with <7g. This crossing point represents the characteristic scale where the power spectrum 
quadrupole changes sign, switching from the linear Kaiser effect to the non-linear velocity 
dispersion effects. In the right panel, the non-linear correction D(k,p) is shown, which gen¬ 
erally increases its departure from unity as ug is increased. Black points are omitted in these 
plots to avoid excessive cluttering, but they are similarly well fitted by the black curves as 
for the colored ones. 

The direct dependence of Pf on the amplitude is better visualized in figure 19, where 
we plot the quantity Pp(/c)fc 3 /( 2n 2 ) for the same models as in figure 18. Analogously to the 
amplitude A 2 {k) in equation (3.5), this is the dimensionless Fourier amplitude of transmission 
fluctuations. Here the black points are also included (with omitted errorbars, which are equal 
to the colored errorbars). At low k, Pf increases very little for p = 0 as og increases, but 
this increase is more pronounced at high p, indicating an increase of [3 with erg. At high k, 
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Figure 19. Amplitude of transmission fluctuations in Fourier space, defined as Ppk 3 / (2n 2 ). The 
power Pp increases slowly with the amplitude of the linear power spectrum parameterized by as, 
except at high k and ^ where it decreases with ag. Black points, for the fiducial model, have their 
errorbars omitted but they are equal to the colored ones. 


Pp continues to increase with the linear amplitude at low /r, but it decreases at high [i. This 
corresponds to the decrease of the wavenumber k na for the crossing of the curves in figure 18 
with increasing a$. 
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Figure 20. Bias and redshift distortion factors as a function of redshift for models SO.76 and SO.64, 
which differ from the fiducial one only for a lower linear amplitude of the mass power spectrum, and 
of the Planck model, which also has a different power spectrum slope. All values are obtained from 
the fits with Di(k,n) fixing <72 = 0. 

The change of bias and redshift distortion factors is shown in figure 20. The very weak 
dependence of Pp on the linear amplitude at low k and low ^1 implies that b T g decreases with 
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the amplitude nearly as fast as cr^ 1 . At high /r the power Pp has a substantial increase with 
us, so the bias factor b Trj decreases more slowly with erg, particularly at high redshift, and as 
a consequence (3 increases with cr 8 , as seen in the left panel. 

6.2 Power spectrum slope 




Figure 21. Comparison of Pf/Pl (left panel) for the fiducial model (with n s = 1) and the Planck 
model (with n s = 0.9624). 

We now compare the model labeled as ’’Planck” in table 3.1, which has parameters 
consistent with the ones measured by the Planck mission, with the fiducial model. Figure 
24 shows Pf/Pl and D(k,fi) for the Planck model as colored curves and points, compared 
to the fiducial model as black curves (with points omitted). The transmission fluctuation 
amplitude, Pp(k)k 3 /(2ir 2 ), is shown also in figure 22, in all cases at z = 2.6. Interpreting the 
differences between these two models is in this case more complicated than in the previous 
subsection, because now both the amplitude and slope of the power spectrum are different. 
In addition the density-temperature relation also differs slightly for these two models as 
specified in table 3.1. 

The fiducial and Planck models have slightly different values of Horn, implying a differ¬ 
ence in their growth factors. It is therefore useful to express the mass fluctuation amplitudes 
on spheres of radius 8 h -1 Mpc at z = 2.6, instead of the present epoch, which have the 
following values: cr 8 (z = 2.6) = (0.3102,0.2910,0.2679) for the fiducial, Planck and SO.76 
models, respectively. The Planck model therefore has a linear amplitude that is close to the 
average of the fiducial and SO.76 models at the scale of <7$. In addition, comparing figures 19 
and 22, we can see that the amplitude of Pp of the Planck model is intermediate between 
that of the fiducial and SO.76 models at the lowest values of k probed by our simulations, 
k ~ 0.2/i/Mpc. This value of k roughly corresponds to the scale at which the parameter 
us measures the amplitude of density fluctuations, and can therefore act as an approximate 
pivot scale for comparing amplitudes of models with different slope. In agreement with this, 
the Planck model has a bias factor b T s intermediate between the values for the fiducial and 
SO.76 models (see figure 20). This bias factor has therefore little sensitivity to the power spec¬ 
trum slope, and depends mostly on the amplitude parameter cr$(z). We caution, however, 
that this bias factor may also depend on the temperature-density relation, and is sensitive to 
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Figure 22. Amplitude of transmission fluctuations for the Planck model (colored curves and points) 
and the fiducial model (black curves and points with errorbars omitted). 


the value of F. The behavior of b TT] and /3 shows a steeper redshift evolution of the Planck 
model compared to the fiducial one with varying amplitude. This different evolution may be 
affected by the lower temperature of the Planck model compared to the fiducial one, as will 
be discussed in §6.4. 

The ratio of the linear power spectra for the Planck model to that of the fiducial model 
varies approximately as A: -0 ' 03 over the range of scales we explore, as caused mostly by the 
difference in the primordial spectral index n s (with a smaller effect arising from the different 
value of Oom h). One might expect a similar ^’-dependence of the ratio of the transmission 
power Pp for the two models, seen in figure 22, but this would generally be misleading. In 
fact, figure 19 shows that the ratio of Pp for the SO.76 and fiducial models has a similar 
fc-dependence at low //, even though these two models have exactly the same Pp{k) except 
for the normalization. This is caused by the increase of the non-linear correction D with as 
at low k and [i. 

In general, the sensitivity of Pp to Pp is relatively weak on scales k > 0.3/i -1 Mpc, 
when non-linear effects start. This means that Pp/Pp tends to behave inversely with Pp, 
as found in §6.1 when analyzing the dependence of Pp on the amplitude, and therefore to 
increase faster with k for the Planck model compared to the fiducial one. Of course, on very 
large scales the constancy of the bias factors implies that Pf is proportional to Pp for a given 
model. At the same time, on very small scales (k > 3/i^ 1 Mpc), the reason for the slower 
decline of Pp/Pp for the Planck model compared to the fiducial one is caused mostly by the 
lower temperature in the Planck simulation (see §6.4). 

Finally, figure 23 shows the non-linear parameters of the function Di(k,fp) obtained 
in our fits for the fiducial, SO.76, SO.64 and Planck models. The parameter k p is mostly 
dependent on the gas temperature and is higher for the Planck model owing to its lower 
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Figure 23. Parameters of the non-linear function Di(fc,/i) for the fiducial model (red squares), 
for the SO.76 and SO.64 models with a decreasing amplitude of the linear power spectrum (green 
triangles and purple rhombi), and for the Planck model, with an intermediate amplitude between the 
fiducial and SO.76 models and a lower n s (blue circles). The sixth panel shows the parameter k na , 
the wavenumber at which curves of different /i cross each other in the diagram of Pp(k, /F)/PL{k), 
computed from equation (4.4). 


temperature, and very similar for the other three models. The parameters q \, a v and k/ v 
mostly follow a progression depending on the amplitude of Pl : as the amplitude increases, 
qi is decreased because the impact of non-linear terms does not increase as fast as A 2 (A;) 
in equation (3.6), and a v and k% v are also reduced. Variations of b v are rather small. The 
sixth panel in the figure shows k na (z), representing the scale at which the curves of Pf/Pl 
at different /r cross each other. As expected, the errors of this parameter are small because 
most of the degeneracy with other parameters is removed. The variation with redshift and 
us confirms what we have described: as the power amplitude increases, k na decreases. 

6.3 Mean transmission fraction and redshift evolution 

The observed evolution with redshift of the non-linear power spectrum is a consequence of 
two different effects: the physical evolution that would be observed at fixed F. and the 
change in Pp that occurs when varying F at a fixed redshift. These two effects are clearly 
separated in figure 24, where the change from z = 2.2 to z = 3 for Pf/Pl is shown for the 
fiducial model at fixed F. In the left panel, we choose the value of F that corresponds to 
z = 2.2, and the standard result for z = 2.2 for the fiducial model is shown as the black 
curves. Colored curves and points show the result at z = 3 when using the value of F that 
corresponds to z = 2.2. Similarly, the black curves in the right panel show our standard 
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Figure 24. Evolution of Pf/Pl with redsliift when the mean transmission is kept fixed at F = 0.851 
(left panel) and F = 0.696 (right panel). The result of the fiducial model is shown in both panels 
at z = 2.2 and z = 3, as black curves (with points omitted) when the true value of F is used, and 
as colored curves and points when F is modified to the value at the other redshift. The observed 
evolution, from the black curves in the left panel to those in the right panel, is a result of both the 
physical redshift evolution at a fixed F, and the variation of F with redshift. 


result at z = 3, while the colored curves and points show the z = 2.2 result for the value 
of F that corresponds to z = 3. The density-temperature relation of our models is nearly 
constant with redshift, so this should have a minimal impact on the variations seen in this 
figure. 

The redshift evolution at fixed F arises basically as a consequence of the varying ampli¬ 
tude of the linear power spectrum. The amplitude crs(z) varies by 4/3.2 = 1.25 from z = 3 
to z = 2.2 (neglecting the influence of the cosmological constant at these high redshifts), so 
the redshift evolution behaves in the same way as the variation with as in figure 18: Pf/Pl 
decreases nearly as fast as ag(z)~ 2 at low k and /i, with an increase of /3 with as(z), and 
the value k na where the quadrupole changes sign decreases (i.e., shifts to larger scales) with 
ag(z). ft is interesting to note that this value of k na is not shifting at all under the large 
change of F from 0.851 in the left panel to 0.696 in the right panel, and depends on the 
redshift only despite the clear variations of the shape of the curves with F. 

The bias and redshift distortion factors are shown in figure 25 for the fiducial model, 
for three different cases. The red squares show their true evolution when F(z) is varied at 
the same time as the redshift. The blue circles show the isolated redshift evolution, when 
keeping F = 0.781 fixed to its value at z = 2.6, and the green triangles show the result when 
varying only F but keeping fixed z = 2.6 (only the redshift z is shown in the horizontal axis 
but the green points are all for z = 2.6, and for F equal to the value at the redshift in the 
axis). 

The redshift evolution of b T s for fixed F is, as mentioned in §6.1, equivalent to the vari¬ 
ation with the power spectrum amplitude, increasing with redshift nearly as b T s oc og^) -1 , 
in agreement with the blue points in the middle panel. When varying F at fixed redshift, 
Pp increases rapidly as F declines but this is mostly due to the relation between bps and 
b T s . The residual variation left for the optical depth bias on F, shown as the green points, 
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Figure 25. Bias and redshift distortion factors for the fiducial model. Red squares show our standard 
result for the redshift evolution, when F(z) varies with redshift as observed. Blue circles show the 
redshift evolution when we fix F = 0.781, equal to the value at z = 2.6. Green triangles show the 
variations at a fixed redshift z = 2.6, varying only F(z) to the value it has at each of our standard 
redshift outputs. Blue circles (green triangles) are shifted slightly right (left) to avoid superposition 
of errorbars. 

is actually of opposite sign: as F decreases, regions of declining overdensity in the IGM are 
responsible for the dominant variations measured in the Lya forest spectra, and they have 
a decreasing physical bias factor. The true evolution we can measure of b T $, shown as the 
red points, is nearly constant with redshift due to a fortuitous cancellation of the effects 
of declining power spectrum amplitude and declining F with redshift. This cancellation is, 
however, not a very precise prediction from our simulations because it is sensitive to the fit 
we use to obtain the bias factors, as can be seen from the redshift evolution of the fiducial 
model bias factors for the free q 2 fit that is shown in figure 14. 

The right panel of figure 25 shows that the usual decrease of b TV with redshift is due to 
the variation of F. In fact, for fixed F. b rrj actually has a slow increase with redshift. The 
reason why the redshift distortion factor (3 generally decreases with redshift is therefore not 
very straightforward. At fixed F, the evolution of (3 has a relatively simple explanation: the 
decreasing power spectrum amplitude with redshift implies an increase of b r s with redshift, 
while b rr] does not change much at fixed F, so (3 decreases with redshift. However, the reason 
why the ratio 6 r); / b T $ has a redshift evolution that is not strongly changed by the variation 
of F is not so simple. 

The non-linear parameters of the fits can be seen in figure 30. Some of the variations 
with redshift correspond to the variations we have seen before with the amplitude in figure 23: 
the parameters q\ , k% v and a v increase with redshift for fixed F. However, these parameters 
also increase with F and the result for the true evolution is a much slower variation. It is 
interesting that the parameter k p also varies with redshift when computed at fixed F and 
that this is also cancelled by the dependence on F, an effect that likely depends on the slope 
of the density-temperature relation. 

6.4 Temperature-Density relation 

To finish this section, we investigate the dependence of Pf/Pl on the temperature-density 
relation. First, we compare in figure 27 two models that have both 7 = 1 (a flatter relation 
than in our fiducial model, which has 7 = 1 . 6 ), when we vary the temperature Tq at the 
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Figure 26. Parameters of the non-linear function Di(fc,/i) for the cases of true redshift evolution 
(red squares), fixing F to its value at z = 2.6 (blue circles), and fixing z = 2.6 but using the value of 
F corresponding to each redshift (green triangles). 


mean density from 10 4 K in model G1T4, shown as colored points and curves, to 10 4 ' 3 K in 
model G1.0, shown as black curves and points. The results are shown at z = 2.6 and the fits 
are for q 2 = 0. A lower temperature reduces the Jeans length and, as expected, the damping 
of the power at small scales is reduced, as seen clearly in figure 27. This change in the small- 
scale power induces a modification of the bias factors on large scales, which as discussed for 
previous cases, depends on the requirement of keeping a fixed value of F. Unfortunately, 
when we compare the difference of the colored and black points at low k in figure 27 with the 
difference between the fits, we see that whereas the points do not indicate any variation of 
f3 with temperature (the difference in the points is constant at low k without an appreciable 
dependence on /x), the fits show a different behavior: they reflect the difference between the 
two models at low n, but at high n the two curves are nearly the same, implying a higher 
j3 for the low-tenrperature models. This results from the need to obtain a good fit to the 
well-determined points at high k , and the large errors at low k, and therefore the implied 
small change of /3 with temperature is not reliable. 

The more reliable modification induced by the gas temperature is the variation of the 
scale k na : a lower gas temperature implies lower velocity dispersions and therefore smaller 
non-linear effects on the power spectrum anisotropy, so the scale of sign reversal of the power 
spectrum quadrupole decreases. Figure 27 shows that k na decreases by ~ 10% when the 
temperature increases by a factor of two. This provides a possible clue to resolve degeneracies 
with (is and F in using the shape of the non-linear power to constrain the IGM temperature 
from future measurements (e.g., by using both the value of k na and the damping rate of the 
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Figure 27. Dependence of Pf/Pl on the gas temperature, at z = 2.6. The two models shown have a 
temperature-density relation slope 7 = 1 , and the gas temperature is twice higher for the black points 
compared to the colored ones. 



k [h/Mpc] 


Figure 28. Dependence of Pf/Pl on the slope of the temperature-density relation at z = 2.6. The 
black points are for the fiducial model, and the colored points show the G1.0 model result, when 
modifying the slope to 7 = 1 (which is the same model as the black curves in figure 27. 
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power at higher k). 

The effect of varying the temperature-density relation slope is shown in figure 28, where 
we compare the G 1.0 model (with 7 = 1 ; colored curves and points) with the fiducial model 
(7 = 1.6; black curves and points), also at 2 = 2.6. The temperature is kept fixed at the mean 
density, at To = 10 4 ' 3 K. There is also a reduction of the damping at small scales in the G1.0 
model, probably because the relevant value of the temperature is at a density above the mean 
when we use F(z = 2.6) = 0.781, where the G1.0 model has the lower temperature. This 
agrees also with the fact that the value of k na decreases with the slope 7 , in the same way as 
it decreases with the temperature as in figure 27. However, the reduction in slope produces 
an increase of the bias factors, contrary to the temperature reduction which decreases them. 
Therefore, sufficiently careful observations of the non-linear shape of the power spectrum 
may help disentangle many of the properties of the IGM. 




Figure 29. Bias and redshift distortion factors for models with varying temperature-density relation, 
T = T 0 (p/p) 7_1 . Three models have a temperature at the mean density To = 10 4 3 K, and slopes 
7 = 1.6 (fiducial, red squares), 7 = 1.3 (G1.3, purple rhombi), and 7=1 (G1.0, blue circles). The 
model G1T4 (green triangles) also has 7=1 and a lower temperature, To = 10 4 K. 

The redshift distortion and bias factors of these models are shown in figure 29. A fourth 
model of intermediate slope between the fiducial and G1.0 models is added. The figure 
confirms the general property we have mentioned above: the density bias factor increases 
with the temperature To and decreases with the slope 7 , but the variation is small. The 
general conclusion that (3 drops with redshift is valid for all models, but the steeper redshift 
evolution of /3 for the low-temperature model G1T4 seen in the left panel is not so reliable 
because of the difference between the fits and the actual simulation points at low k described 
in figure 27. This detailed question can only be resolved by examining more simulations of 
different temperature-density relations on large boxes. 

The non-linear parameters obtained from the same fits are shown in figure 30. The 
parameter k p increases as expected when To is reduced for the G1T4 model (green triangles), 
corresponding to a decreasing Jeans scale, and the flattening of the slope also increases k p 
even for the models of fixed To, as mentioned above. However, this variation of k p with To 
has a large dependence on redshift. A much clearer effect of reducing To is the reduction 
of ky v and a v , which is not affected by the slope 7 . The model dependence of k na , which 
has much smaller errors due to the much smaller degeneracy with other parameters, shows 
clearly the effects that we have referred to earlier: the wavenumber k na decreases with the 


- 44 - 



































Figure 30. Parameters of the non-linear function Z?i(fc,//) for the fits to the same four models with 
different temperature-density relation as in figure 29. 


temperature To, and decreases also to a smaller extent when the slope 7 is increased. There 
are also interesting changes of the power-law index b v for the //-dependence with the slope 
7, suggesting again that a precise measurement of the multiple characteristics of the non¬ 
linear power spectrum in future observational studies can provide a precise diagnostic of the 
evolution of the IGM. 

7 Discussion 

The study of the Lya forest has reached an age of maturity. The observations are providing 
precision measurements of the transmission correlation function or power spectrum: the 
one-dimensional power spectrum is measured to very high accuracy [46], and the three- 
dimensional correlation function has been measured on large-scales with a new method that 
treats the distortion introduced by quasar continuum fitting [3]. In this context, reliable 
theoretical predictions are needed to confront the large amount of observational results that 
are becoming available from large surveys [34]. 

In this work, we have introduced a new formula for fitting the numerical results on the 
Lya power spectrum from simulations of structure formation, designated as D\ in equation 
(3.6). This formula was obtained from a guess for how the non-linear correction should 
depend on the fluctuation amplitude in Fourier space, A 2 (k), but it is not derived from any 
perturbation theory. It is based on trial-and-error experimentation to see which formula 
best fits the numerical results, and on the previous work by M03, where a //-dependence of 
the form D oc exp[— A(k)fi bv ] was proposed. We have found in this work that in fact, this 
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/i-dependence of the non-linear correction is a surprisingly good fit to the simulation results 
(figure 6 ), but we have no analytic explanation for this. In the future, it would be desirable to 
attempt developing a second-order perturbation theory for modeling the non-linear coupling 
of a given Fourier mode to all other modes, which suggests that D(k, fi) should be given by 
an integral over all possible wavevectors k 7 of the amplitude product A(k , )A(k — k') times 
a coupling strength [37], instead of our simple term proportional to A 2 (k). 

The form of the non-linear power spectrum is clearly predicted by our simulations, and 
basically agrees with that found in M03. We have characterized it in terms of the wavenum¬ 
ber, k na , at which non-linear anisotropy due to velocity dispersion and thermal broadening 
starts dominating over the linear anisotropy of large-scale velocity flows, and the curves of the 
ratio Pp(k, ii)/Pi,{k) at different /i cross each other. At higher wavenumbers the power falls 
fast for high [i due to thermal broadening, and more slowly at low /j due to the Jeans scale 
of gas at the densities of typical absorption systems, with a characteristic smoothing scale 
given by our k p parameter. As the precision and reliability of the hydrodynamic simulations 
improves, more accurate predictions on the precise shape of the non-linear power spectrum 
should be possible. This requires only analyzing enough simulations with large boxes and 
sufficient resolution to reduce the systematic effects we discuss in §5. The details of the shape 
of Pp(A:, /i) should probably allow to measure independently the amplitude and slope of the 
power spectrum PrXk) and the temperature-density relation, as described in § 6 . 

The numerical results also make a prediction of the two bias factors of the Lya forest. 
The form of the linear power spectrum is that in equation (2.5), but the values of the two 
bias factors depend on all the model parameters that affect the non-linear scales. We have 
introduced the bias factors based on the effective optical depth in equation (2.4), which are 
physically interpretable in the same way as the bias of galaxy populations. Our predictions for 
these bias factors are not highly accurate because of the sampling variance of our simulations 
with limited box size and, to a smaller extent, the limited resolution. Moreover, their fitted 
values depend on the formula that is used to fit the whole shape of the non-linear power 
spectrum, as illustrated by the different results we obtain when the t /2 parameter is left 
free or is set to zero. An improvement over the work presented here should be to consider 
also the cross-power spectrum of 5f with the initial conditions of the simulation, which can 
give an improved estimate of the bias factors. Nevertheless, as explained in §5 we believe 
that our estimates of the bias factors are reliable within ~ 10%. The basic result is that 
the density bias factor is ~ 0.6 for the standard CDMA model, decreases with the power 
spectrum amplitude as cr^ 1 , and is close to constant with redshift as a result of the opposite 
effects of the decreasing power spectrum amplitude and decreasing F with redshift. The 
peculiar velocity gradient bias, which is unity for an isotropic tracer, is predicted to be close 
to but below unity, and decreasing with redshift. The resulting redshift distortion factor, /?, 
is ~ 1.4 at z = 2.2, decreases with redshift and increases with the amplitude <rs. 

The difference of b rri from unity is an important prediction of the hydrodynamic simu¬ 
lations we are studying. If the absorption line features in the Lya forest were arising from a 
series of discrete clouds that are small compared to the scales where cosmological structures 
are entering the non-linear regime, then a change in the large-scale peculiar velocity gradient 
r/ (see equation 2 . 2 ) should simply compress or stretch the fixed absorption profiles into a 
smaller or larger redshift range, depending on the sign of the change in rj. Even though 
the lines increase their blending when tj increases, the effective optical depth still varies as 
(1 — r/) _1 , and so b T1j must be exactly unity. 

It is only when the absorbers have internal velocity flows that are systematically aligned 
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with the large-scale peculiar velocity gradients that they can respond differently to a varying 
7 / at fixed 5. This is the situation when the absorbers in the Lyct forest correspond to 
collapsing structures in the cosmic web of the IGM: if 7 / increases, the absorption lines in 
the forest are compressed into a narrower redshift range, but they also become internally 
narrower and therefore more saturated owing to their own internal velocity gradients, which 
are correlated with those on large-scales. As a consequence, the effective optical depth they 
produce increases slower than 77 , and b TT1 is less than one. This difference from unity should 
be more pronounced when the absorbers are low-density structures following the velocity 
gradients on large-scales, and less pronounced when the dominant absorbing structures are 
due to high-density gas in more virialized regions, where the velocity gradients and dispersions 
are randomized and less correlated with the surrounding large-scale structure. As F increases 
closer to unity (i.e., there is less mean absorption), the dominant absorption lines correspond 
to more overdense gas, and b rr) should be closer to one. This agrees with the behavior we 
find in our simulations (see green points in right panel of figure 25). We also find the value of 
b Trj is substantially larger than 6 r p in equation (2.8), as shown in the right panel of figure 4, 
which would be the value predicted if the internal peculiar velocities of the absorbers could 
be modeled linearly [51]. 

We now compare our results on the non-linear power spectrum with the ones obtained 
by M03. The general shape of Pf/Pl is in very good agreement, as can be seen by comparing 
Figure 9 in M03 with several of our figures, for example, figure 18. For a more quantitative 
comparison, it is useful to compute first the power spectrum normalization at z = 2.25, the 
only redshift at which results were shown in M03. Our Planck model, which has a very 
similar spectral index and gas temperature as the M03 model, has o&{z = 0) = 0.834 and 

= 0.3175, implying ag(z = 2.25) = 0.321, and crg(z = 2.6) = 0.291. The model used in 
M03 had ag(z = 0) = 0.79 and Q m o = 0.4, implying ag(z = 2.25) = 0.290, very close to the 
amplitude of our fiducial model at z = 2.6. Figure 9 of M03 also shows the feature that the 
curves at different /i cross each other nearly at the same point, at k na ~ 4/i/Mpc. This is 
in fact nearly the same value found for our Planck model at z = 2.6, as seen in figure 24, in 
agreement with our finding that this scale depends mostly on the power spectrum amplitude 
(there is also a dependence on gas temperature, which is similar in our Planck model and the 
one in M03). The way the power drops at k > k na for the two models is also fairly similar. 

At low k, the value found for the redshift distortion parameter in M03 of (5 = 1.58 is 
higher than ours. Taking into account the lower power spectrum normalization and lower 
value of F(z = 2.25) = 0.8 used in M03 compared to our fiducial model, our prediction for the 
model used in M03 would be /3 ~ 1.32 when using our q -2 = 0 fits, as derived by interpolating 
results in our figures 20 and 25. Using our free (72 fits generally reduces the predicted value 
of (5 by a further few percent. There is therefore a substantial discrepancy in the redshift 
distortion factor. The density bias factor given in Table 1 of M03 of = 0.0173 is, on 
the other hand, in very good agreement with our 52 = 0 fits: the implied b T s = 0.59 for the 
mean transmission F = 0.8 used in M03 is nearly equal to our predicted value for the Planck 
model (see figure 20), which has the same amplitude at z = 2.6 and a similar value of F at 
this redshift. This suggests that the discrepancy in (3 is related to a higher prediction for b rri 
from M03 compared to our models. However, if we use our free q 2 fits, then our value of b T s 
increases by ~ 12 %, and then the discrepancy in f3 is mostly due to the low prediction for 
b T s from M03. 

The reason for this discrepancy and the correct prediction for /3 and b T s will only be 
resolved with the analysis of more simulations on large boxes. Understanding the systematic 
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errors and the accuracy of the prediction of b TT} from hydrodynamic simulations is especially 
important because if this value can be predicted to an accuracy better than 1%, then a 
precise measurement of (3 and b T s yields an observational determination of the growth factor 
logarithmic derivative, /(H), from equation (2.6), and therefore a fundamental test of the 
presence and evolution of dark energy in the Universe at high redshift. 

We also find that the variations of the bias factors with the model parameters specified 
in Table 1 of M03 are in broad agreement with our results described in §6. 

7.1 Comparison to observations 

There are several observations of the Lya forest power spectrum that our predictions can be 
compared to. First, on the linear regime, the measurements of the Lya autocorrelation on 
large scales yielded a first detection of redshift distortions in [52] using BOSS. More recently, 
an improved method that effectively corrects for distortions introduced by continuum fitting 
has been presented in [3], allowing for accurate measurements of the bias factors. Apart from 
this, the one-dimensional power spectrum is a projection of the full redshift space transmission 
power spectrum which includes non-linear effects, and was measured first in [41], and then in 
[46] with the higher accuracy allowed by BOSS. The full redshift-space power spectrum down 
to small scales, where non-linear effects are important, can also be measured from the BOSS 
data and from close pairs of quasars where correlations can be measured down to the smallest 
scales. Results have been presented constraining the Jeans scale [33], but a complete analysis 
of the available data in terms of the models for the full shape of the three-dimensional power 
spectrum is yet to be done. 

Here, we make only a brief comparison with the measured bias factors in [3], leaving 
the comparison with the one-dimensional results and other measurements of the non-linear 
regime for future papers. The result of [3], obtained at a mean redshift z = 2.3 from 
a fit restricted to large scales (from 40/i _1 Mpc to 160/i _1 Mpc) is f3 = 1.39 ± 0.11, and 
6^(1 + j3) = —0.374 ± 0.007; the second quantity is chosen because it has the smallest 
marginalized observational error. 

The measured value of /3 is in excellent agreement with our predictions. Our Planck 
model is the one that uses parameters consistent with present constraints from the Cosmic 
Background Radiation, and predicts j3 = 1.4 at z = 2.3 for the Q 2 = 0 fit, which drops to 
1.3 for free q 2 ■ We note from figures 12 and 14, however, that this prediction is subject 
to a substantial systematic error due to box size and resolution effects. This value also 
has some sensitivity to other poorly determined parameters of the IGM, in particular to F, 
but we see in figure 25 (green triangles) that we expect a weak variation of /3 in the range 
0.85 < F < 0.78. 

Comparing the predicted and measured values of b T s is more model-dependent. The 
results of [3] are obtained using their fiducial model, with H m o = 0.27 and a§(z = 0) = 0.79, 
implying ag(z = 2.3) = 0.307. While the measurement of /3 is not affected by this fiducial 
model used to fit the data, the density bias factor is affected because the Lya transmission 
power that is actually observed is proportional to [as(z)bFs] 2 ■ Our Planck model has a 
slightly higher amplitude, crs(z = 2.3) = 0.313. Correcting for this, if our Planck model had 
been used for the data analysis, the [3] result would then be bps(3 + /3) = —0.367 ± 0.007. 
The theoretical prediction of our Planck model at z = 2.3 using our q 2 = 0 fit, and our 
standard value of F(z = 2.3) = 0.836, is b T , 5 = 0.586 and f3 = 1.40, implying bps = —0.105, 
and bFs( 1 + /3) = —0.253. This is substantially smaller than the measured value. 
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There are two systematic errors that dominate the uncertainty in this comparison. The 
first is the error of our theoretical predictions, dominated by the uncertainty in fitting the 
bias parameters from the low-fc modes. If we use the fit that includes the q 2 free parameter 
in equation (3.6), the predicted density bias factor increases to b T s = 0.68 (see Table 4 
in Appendix B), and with (3 = 1.3, bp$( 1 + ,8) = —0.28. This is still substantially below 
the observed value. The second systematic error is related to the value of F, which is 
poorly determined by observations. Assuming that F = 0.8 at z = 2.3 (instead of the value 
F = 0.836 derived from equation 3.1), we find from the results in figure 25 that the predicted 
bias for Q 2 = 0 would change to b T s — 0.56, implying bps = —0.125 and bps( 1 + f3) = —0.300. 
If we combine both systematic errors (assuming the higher b T $ from our free q 2 fits and a low 
mean transmission F = 0.8 at z = 2.3) yields a value bps( 1 + /3) = —0.334, closer but still 
not consistent with the observational determination of [3]. 

It remains to be seen if the mean transmission from the Lya forest is in fact substan¬ 
tially lower than was found by [31], or if there is a large discrepancy between the predicted 
value of b T s from the simulations analyzed in this paper and the observations. The recent 
determinations of F by [1] give a value similar to that of [31]. It is also possible that the 
value of the bias factors is altered by the presence of intensity fluctuations in the ionizing 
background. The need to understand the observed values of the bias factors highlights the 
importance of obtaining more precise predictions from cosmological simulations. 

8 Conclusions 

The analysis of a variety of simulations of the Lya forest we have presented highlights the fact 
that the detailed form of the Lya transmission power spectrum in redshift space is in principle 
predictable from a well-defined theory of the IGM. The theory assumes that photoionization 
by the cosmic ultraviolet background and the gravitational evolution of structure that arises 
from primordial fluctuations in the Cold Dark Matter scenario determine the properties of 
the IGM and the Lya power spectrum. Detailed observations in the future should provide 
tests of the validity of this simple theory, and investigate to what extent the effects of galactic 
winds, quasar jets, fluctuations of the ionizing background intensity, and inhomogeneities of 
the density-temperature relation arising from reionization, have an impact on the observable 
power spectrum and other characteristics of the Lya forest. 

Much work remains to be done both on the theoretical and data analysis front to fully 
confront the model predictions with observations. The existing data from BOSS and close 
quasar pairs that have been individually observed should provide powerful measurements 
of the full Lya power spectrum, beyond the constraints obtained from the one-dimensional 
analysis [46, 47]. This will also be complemented by new surveys of absorption spectra 
that will improve on the BOSS results. At the same time, the theory needs to be further 
developed with the analysis of more simulations on large boxes, as they are made possible by 
modern computational technology. An improved understanding of the effects of variations 
on the physical model of the IGM is necessary before one can use the non-linear Lya power 
spectrum to constrain fundamental parameters in cosmology related, for example, to neutrino 
masses or other modifications of the dark matter sector. 

To make the results of our simulations accessible for comparing to future simulations or 
observations, we are providing in Appendix B the fit parameters to the most useful models 
that are presented in this paper. At the same time, for more direct comparisons to the results 
of our simulations, we are also providing the results for Pp(k,/j of all our models, with the 
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bins that we have used to obtain all our fits. These data are publicly available at GibHub 3 
repositories (https : //github. com/andreuandreu/3D_Power_spectrum_modes_tables). 
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A Dependence of the Bias Factors on the Minimum Error parameter 

The fits to the non-linear Lya power spectrum we have obtained depend on several technical 
parameters for binning the power spectrum and evaluating the y 2 function, which are de¬ 
scribed in §3.3. We have checked that our results for the best fits are not strongly dependent 
on the way the Fourier modes are binned. As long as our parameter kt is not too small, 
the rebinning of Fourier modes should not modify the fit. However, the e parameter that is 
inserted in the expression for the error in equation 3.2 inevitably affects our fits. The intro¬ 
duction of this parameter is necessary to reduce the weight of the high-A; modes for our fits, 
because the number of modes below k increases as k 3 . We note that if the total number of 
bins nfe at k > kt to evaluate the power Pp is modified (in our fits this was fixed to n& = 256), 
then our fits are modified depending on the paremeter en h 3//2 , because the fits depend on 
the effective weight assigned to modes in each fixed region of Fourier space. 

The choice of the parameter e cannot be made in a very objective way. As e is increased, 
an increasing discrepancy in the high-A; modes is allowed, taking into account that our fit¬ 
ting formula is not a perfect theory that should precisely match the results obtained from 
numerical simulations, which have very small statistical errors at high k that are at some 
point smaller than the systematic errors. Reducing the weight of the high-A; modes allows 
for a better fit of the bias factors from the low-A; modes. 




e e 

Figure 31. Dependence of the redshift distortion factor /? and the density bias factor b T $ on the 
parameter e determining a minimum error in our fit with the D\ formula with <72 = 0. Blue squares 
(red circles) show the result for our fiducial (L120) model, at z = 2.2. 

Therefore, we look into the variation of the redshift distortion factor and the density 
bias factor in the two panels of figure 31. As e is increased, the value of f3 converges just 
near the value we have chosen in this paper e = 0.05. The value of b T s continues to increase 
above this value of e, by ~ 10% as e is increased to 1. Interestingly, the bias factor reaches a 
value b T s — 0.61 at very high e, very close to the result we obtain in our fits with free q -2 for 
these models (see figure 14). This suggests that our models with free q 2 yield more accurate 
values of the bias factors, that are less affected by the need to correctly fit the numerous data 
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points for Pp at high k from our simulations. Nevertheless, this will need to be tested with 
results from more simulations on large boxes than are analyzed in this paper. 
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B Tables 


Table 2. Bias parameters for simulations with different resolution and box size, for the <72 = 0 fit. 


z 

P 

^tS 

b T r] 

Fiducial 

3.0 

1.205 ± 0.049 

0.5546 ± 0.0086 

0.681 ±0.028 

2.8 

1.284 ± 0.052 

0.5577 ± 0.0083 

0.733 ± 0.030 

2.6 

1.343 ± 0.055 

0.5588 ± 0.0084 

0.771 ± 0.032 

2.4 

1.385 ± 0.056 

0.5574 ± 0.0080 

0.796 ± 0.032 

2.2 

1.405 ± 0.061 

0.5536 ± 0.0079 

0.807 ± 0.034 

P1024 

3.0 

1.205 ± 0.050 

0.5566 ± 0.0979 

0.684 ±0.119 

2.8 

1.283 ± 0.052 

0.5597 ± 0.0084 

0.735 ± 0.030 

2.6 

1.340 ± 0.056 

0.5611 ±0.0084 

0.772 ± 0.033 

2.4 

1.381 ± 0.059 

0.5599 ± 0.0082 

0.797 ±0.034 

2.2 

1.403 ± 0.061 

0.5556 ± 0.0080 

0.809 ± 0.035 

L120 

3.0 

1.195 ±0.039 

0.5710 ±0.0057 

0.695 ± 0.023 

2.8 

1.283 ± 0.040 

0.5768 ± 0.0058 

0.757 ±0.024 

2.6 

1.355 ± 0.042 

0.5808 ± 0.0057 

0.808 ± 0.025 

2.4 

1.411 ± 0.042 

0.5820 ± 0.0059 

0.847 ±0.025 

2.2 

1.443 ± 0.044 

0.5789 ± 0.0051 

0.867 ±0.025 

R384C 

3.0 

1.210 ±0.060 

0.5661 ± 0.0094 

0.698 ± 0.035 

2.8 

1.290 ± 0.060 

0.5714 ± 0.0097 

0.754 ±0.036 

2.6 

1.360 ± 0.060 

0.5732 ± 0.0086 

0.800 ± 0.035 

2.4 

1.420 ± 0.050 

0.5715 ± 0.0065 

0.837 ±0.029 

2.2 

1.455 ± 0.064 

0.5588 ± 0.0085 

0.844 ±0.037 

R384 

3.0 

1.180 ±0.060 

0.5708 ± 0.0100 

0.687 ±0.036 

2.8 

1.270 ± 0.060 

0.5743 ± 0.0096 

0.746 ± 0.036 

2.6 

1.330 ± 0.060 

0.5760 ± 0.0086 

0.787 ±0.035 

2.4 

1.380 ± 0.070 

0.5758 ± 0.0087 

0.820 ± 0.041 

2.2 

1.420 ± 0.070 

0.5703 ± 0.0099 

0.841 ± 0.041 

R640 

3.0 

1.251 ±0.052 

0.5409 ± 0.0953 

0.690 ±0.120 

2.8 

1.316 ± 0.054 

0.5456 ± 0.0085 

0.735 ±0.031 

2.6 

1.364 ± 0.054 

0.5482 ± 0.0082 

0.768 ±0.031 

2.4 

1.397 ±0.056 

0.5487 ± 0.0079 

0.791 ±0.032 

2.2 

1.421 ±0.061 

0.5461 ± 0.0085 

0.805 ± 0.034 
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Table 3. Non-linear fit parameters for the same simulations as in Table 2. 


z 

Qi 

k p [/i/Mpc] 

k%° [h/Mpc] a - 

CLy 

b v 

h 

^na 

Fiducial 

3.0 

0.648 ± 0.042 

13.1 ±0.4 

1.007 ±0.103 

0.627 ±0.033 

1.66 ±0.01 

4.09 ± 0.04 

2.8 

0.644 ± 0.037 

13.4 ±0.5 

0.979 ± 0.091 

0.610 ±0.030 

1.64 ± 0.01 

3.89 ± 0.03 

2.6 

0.652 ± 0.034 

13.6 ±0.5 

0.970 ± 0.084 

0.590 ± 0.028 

1.61 ± 0.01 

3.65 ± 0.04 

2.4 

0.666 ± 0.030 

13.5 ±0.5 

0.963 ± 0.076 

0.561 ± 0.027 

1.58 ±0.01 

3.39 ±0.03 

2.2 

0.677 ±0.027 

13.3 ±0.5 

0.961 ± 0.070 

0.533 ± 0.025 

1.54 ± 0.01 

3.11 ±0.03 

P1024 

3.0 

0.644 ± 0.042 

13.2 ±0.4 

0.997 ±0.100 

0.634 ± 0.032 

1.66 ±0.02 

4.09 ± 0.04 

2.8 

0.642 ± 0.037 

13.5 ±0.6 

0.973 ± 0.091 

0.617 ±0.030 

1.63 ±0.02 

3.89 ± 0.03 

2.6 

0.650 ± 0.035 

13.7 ±0.5 

0.964 ± 0.084 

0.597 ±0.029 

1.60 ±0.02 

3.65 ± 0.04 

2.4 

0.665 ± 0.031 

13.6 ±0.5 

0.959 ± 0.075 

0.567 ±0.026 

1.57 ±0.02 

3.39 ±0.03 

2.2 

0.679 ± 0.028 

13.4 ±0.5 

0.956 ± 0.069 

0.537 ±0.025 

1.53 ±0.01 

3.11 ±0.03 

L120 

3.0 

0.661 ± 0.032 

11.3 ±0.2 

0.803 ± 0.064 

0.460 ± 0.026 

1.55 ±0.01 

4.01 ±0.03 

2.8 

0.635 ± 0.028 

12.0 ±0.3 

0.782 ± 0.057 

0.471 ± 0.024 

1.56 ±0.01 

3.80 ± 0.03 

2.6 

0.626 ± 0.026 

12.7 ±0.4 

0.772 ± 0.053 

0.469 ± 0.023 

1.56 ±0.01 

3.56 ± 0.02 

2.4 

0.633 ± 0.023 

13.2 ±0.4 

0.771 ± 0.048 

0.456 ± 0.021 

1.54 ± 0.01 

3.29 ± 0.02 

2.2 

0.643 ± 0.020 

13.5 ±0.6 

0.778 ± 0.046 

0.439 ± 0.020 

1.52 ±0.01 

3.00 ± 0.02 

R384C 

3.0 

0.592 ±0.051 

14.0 ±0.6 

0.786 ±0.104 

0.504 ± 0.042 

1.57 ±0.01 

4.15 ±0.05 

2.8 

0.578 ± 0.045 

14.9 ±0.8 

0.774 ± 0.095 

0.509 ± 0.038 

1.58 ±0.01 

3.92 ± 0.05 

2.6 

0.584 ±0.039 

15.5 ±0.9 

0.773 ± 0.083 

0.500 ± 0.036 

1.58 ±0.02 

3.67 ±0.04 

2.4 

0.602 ± 0.032 

15.9 ± 1.1 

0.778 ± 0.064 

0.478 ± 0.028 

1.55 ±0.01 

3.40 ± 0.03 

2.2 

0.644 ± 0.027 

15.4 ±0.8 

0.822 ± 0.060 

0.458 ± 0.024 

1.52 ±0.02 

3.12 ±0.03 

R384 

3.0 

0.586 ±0.051 

12.9 ±0.5 

0.921 ± 0.120 

0.617 ±0.041 

1.66 ±0.02 

4.13 ±0.05 

2.8 

0.584 ±0.043 

13.3 ±0.5 

0.894 ± 0.104 

0.600 ± 0.038 

1.64 ± 0.02 

3.93 ± 0.05 

2.6 

0.594 ± 0.041 

13.6 ±0.6 

0.885 ± 0.097 

0.578 ± 0.037 

1.61 ± 0.01 

3.68 ± 0.04 

2.4 

0.612 ± 0.038 

13.7 ±0.7 

0.877 ±0.089 

0.547 ±0.034 

1.58 ±0.01 

3.42 ± 0.03 

2.2 

0.634 ± 0.030 

13.5 ±0.6 

0.880 ± 0.079 

0.514 ±0.032 

1.54 ± 0.01 

3.13 ±0.03 

R640 

3.0 

0.642 ± 0.042 

14.9 ±0.8 

1.023 ±0.111 

0.651 ±0.034 

1.71 ±0.03 

4.14 ±0.04 

2.8 

0.633 ± 0.038 

15.1 ±0.7 

0.998 ± 0.094 

0.636 ± 0.030 

1.68 ±0.02 

3.95 ± 0.04 

2.6 

0.636 ± 0.034 

15.1 ±0.7 

0.984 ±0.085 

0.614 ±0.028 

1.65 ±0.02 

3.71 ± 0.04 

2.4 

0.645 ± 0.030 

14.8 ±0.6 

0.969 ± 0.077 

0.582 ± 0.027 

1.61 ±0.01 

3.45 ± 0.02 

2.2 

0.656 ± 0.027 

14.4 ±0.7 

0.950 ±0.071 

0.543 ± 0.026 

1.57 ±0.02 

3.16 ±0.03 
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Table 4. Bias parameters for D\ with q 2 set free. 


z 

P 

brS 

b T r] 

Fiducial 

3.0 

1.179 ±0.061 

0.6010 ± 0.0126 

0.722 ± 0.039 

2.8 

1.256 ± 0.049 

0.6096 ± 0.0132 

0.784 ±0.033 

2.6 

1.292 ± 0.069 

0.6206 ± 0.0146 

0.823 ± 0.046 

2.4 

1.321 ± 0.067 

0.6286 ± 0.0130 

0.856 ± 0.044 

2.2 

1.347 ±0.047 

0.6306 ± 0.0099 

0.882 ±0.031 

R384C 

3.0 

1.200 ± 0.067 

0.5981 ± 0.0194 

0.732 ± 0.045 

2.8 

1.288 ± 0.068 

0.6043 ± 0.0137 

0.797 ± 0.044 

2.6 

1.355 ± 0.070 

0.6096 ±0.0116 

0.848 ± 0.044 

2.4 

1.399 ± 0.078 

0.6147 ±0.0144 

0.887 ±0.050 

2.2 

1.426 ± 0.074 

0.6144 ± 0.0126 

0.909 ± 0.047 

L120 

3.0 

1.204 ± 0.032 

0.5808 ± 0.0098 

0.713 ±0.021 

2.8 

1.296 ± 0.028 

0.5962 ± 0.0104 

0.791 ±0.021 

2.6 

1.365 ± 0.053 

0.6092 ± 0.0131 

0.854 ±0.036 

2.4 

1.409 ± 0.052 

0.6202 ± 0.0100 

0.901 ± 0.034 

2.2 

1.431 ± 0.036 

0.6249 ± 0.0105 

0.928 ± 0.026 

L80 

3.0 

1.120 ±0.040 

0.6169 ±0.0185 

0.704 ±0.032 

2.8 

1.210 ±0.030 

0.6251 ±0.0115 

0.774 ± 0.023 

2.6 

1.280 ± 0.030 

0.6344 ±0.0117 

0.834 ±0.024 

2.4 

1.340 ± 0.030 

0.6375 ± 0.0078 

0.881 ±0.021 

2.2 

1.350 ± 0.030 

0.6474 ± 0.0120 

0.907 ±0.024 

Planck 

3.0 

1.072 ± 0.059 

0.6462 ± 0.0185 

0.657 ±0.045 

2.8 

1.155 ±0.065 

0.6573 ± 0.0167 

0.723 ± 0.049 

2.6 

1.254 ± 0.044 

0.6601 ± 0.0109 

0.791 ± 0.034 

2.4 

1.290 ± 0.074 

0.6784 ± 0.0154 

0.840 ± 0.056 

2.2 

1.316 ±0.057 

0.6751 ± 0.0141 

0.858 ± 0.044 
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Table 5. Non-linear fit parameters for D\ with q 2 set free. 


z 

Qi 

Q2 

k p [h/Alpc} 

k% v [h/Mpc] av 

d v 

by 

Fiducial 

3.0 

0.104 ±0.121 

0.444 ± 0.093 

10.1 ±0.5 

0.516 ±0.114 

0.248 ± 0.072 

1.66 ±0.02 

2.8 

0.086 ±0.110 

0.417 ±0.078 

9.9 ±0.5 

0.493 ± 0.087 

0.217 ±0.069 

1.63 ±0.02 

2.6 

0.068 ±0.101 

0.390 ± 0.063 

9.6 ±0.5 

0.483 ± 0.052 

0.190 ±0.056 

1.61 ±0.02 

2.4 

0.057 ±0.085 

0.368 ± 0.049 

9.2 ±0.3 

0.480 ± 0.028 

0.156 ±0.045 

1.57 ±0.02 

2.2 

0.090 ±0.052 

0.316 ±0.027 

8.9 ±0.3 

0.493 ±0.016 

0.145 ±0.031 

1.54 ±0.01 

R384C 

3.0 

0.207 ±0.208 

0.320 ±0.159 

11.2 ± 1.0 

0.468 ±0.165 

0.208 ±0.120 

1.57 ±0.01 

2.8 

0.202 ±0.132 

0.289 ± 0.094 

11.4 ±0.8 

0.462 ±0.111 

0.207 ±0.093 

1.57 ±0.01 

2.6 

0.203 ± 0.073 

0.267 ±0.041 

11.2 ±0.4 

0.467 ±0.035 

0.198 ±0.042 

1.57 ±0.01 

2.4 

0.208 ± 0.083 

0.246 ± 0.035 

11.0 ±0.3 

0.475 ± 0.019 

0.182 ±0.027 

1.54 ±0.01 

2.2 

0.233 ± 0.062 

0.213 ±0.023 

10.8 ±0.2 

0.499 ±0.015 

0.181 ± 0.019 

1.51 ±0.01 

L120 

3.0 

0.529 ±0.115 

0.117 ±0.102 

10.6 ±0.6 

0.684 ±0.196 

0.362 ± 0.085 

1.54 ±0.01 

2.8 

0.398 ±0.112 

0.190 ±0.087 

10.5 ±0.6 

0.583 ±0.167 

0.290 ± 0.076 

1.55 ±0.01 

2.6 

0.316 ±0.105 

0.226 ± 0.070 

10.3 ±0.7 

0.529 ±0.118 

0.233 ± 0.074 

1.54 ±0.01 

2.4 

0.264 ±0.091 

0.241 ± 0.060 

9.9 ±0.6 

0.501 ± 0.069 

0.183 ±0.065 

1.53 ±0.01 

2.2 

0.245 ± 0.086 

0.232 ±0.051 

9.6 ±0.6 

0.497 ± 0.046 

0.153 ±0.059 

1.51 ±0.01 

L80 

3.0 

0.144 ±0.232 

0.430 ± 0.202 

8.9 ±0.8 

0.514 ± 0.286 

0.212 ±0.158 

1.61 ±0.01 

2.8 

0.130 ±0.133 

0.394 ±0.108 

9.0 ±0.6 

0.490 ±0.108 

0.190 ±0.093 

1.59 ±0.01 

2.6 

0.117 ±0.110 

0.367 ± 0.080 

8.9 ±0.5 

0.480 ± 0.059 

0.167 ±0.071 

1.56 ±0.01 

2.4 

0.147 ±0.082 

0.320 ± 0.054 

8.8 ±0.4 

0.493 ± 0.041 

0.161 ± 0.054 

1.53 ±0.01 

2.2 

0.120 ±0.083 

0.308 ± 0.047 

8.5 ±0.4 

0.484 ± 0.021 

0.123 ±0.050 

1.50 ±0.01 

Planck 

3.0 

0.0004 ± 0.2241 

0.787 ±0.220 

11.8 ± 1.2 

0.552 ±0.201 

0.178 ±0.108 

1.64 ±0.03 

2.8 

-0.0168 ±0.1686 

0.711 ±0.143 

12.0 ± 1.0 

0.548 ±0.126 

0.203 ± 0.073 

1.65 ±0.03 

2.6 

0.0004 ±0.1241 

0.651 ±0.113 

11.8 ± 1.0 

0.545 ±0.115 

0.195 ±0.073 

1.64 ±0.02 

2.4 

-0.0020 ±0.1108 

0.623 ± 0.069 

10.9 ±0.5 

0.517 ±0.027 

0.152 ±0.035 

1.62 ±0.02 

2.2 

0.0405 ± 0.0890 

0.534 ± 0.049 

10.6 ±0.3 

0.534 ±0.022 

0.139 ±0.029 

1.60 ±0.01 
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Table 6. Bias parameters for Eulerian, Lagrangian and Planck simulations, for <72 = 0. 


z 

P 

^tS 

b T r) 

Fiducial 

3.0 

1.21 ±0.049 

0.555 ± 0.0086 

0.681 ±0.028 

2.8 

1.28 ±0.052 

0.558 ± 0.0083 

0.733 ± 0.030 

2.6 

1.34 ± 0.055 

0.559 ± 0.0084 

0.771 ± 0.032 

2.4 

1.39 ±0.056 

0.557 ±0.0080 

0.796 ± 0.032 

2.2 

1.40 ±0.061 

0.554 ± 0.0079 

0.807 ± 0.034 

Eulerian 

2.9 

1.29 ±0.073 

0.557 ±0.011 

0.737 ±0.042 

2.6 

1.51 ±0.084 

0.561 ±0.010 

0.870 ± 0.048 

2.3 

1.63 ±0.083 

0.557 ±0.011 

0.941 ± 0.047 

Lagrangian 

3.0 

1.03 ±0.055 

0.555 ±0.011 

0.583 ± 0.032 

2.8 

1.14 ±0.059 

0.561 ±0.011 

0.653 ± 0.034 

2.6 

1.24 ± 0.064 

0.563 ±0.011 

0.720 ± 0.037 

2.4 

1.33 ±0.069 

0.558 ±0.010 

0.769 ± 0.039 

2.2 

1.47 ±0.070 

0.556 ± 0.0097 

0.854 ± 0.04 

Planck 

3.0 

1.12 ±0.056 

0.588 ± 0.0096 

0.627 ±0.036 

2.8 

1.22 ±0.055 

0.592 ± 0.0094 

0.689 ± 0.036 

2.6 

1.30 ±0.057 

0.591 ± 0.0093 

0.736 ± 0.037 

2.4 

1.39 ±0.058 

0.591 ± 0.0095 

0.788 ± 0.038 

2.2 

1.42 ± 0.058 

0.581 ± 0.0092 

0.796 ± 0.038 
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Table 7. Non-linear fit parameters for the Eulerian, Lagrangian and Planck simulation, for <72 = 0. 


z 

Qi 

k p [h/Mpc] 

ht? [h/Mpc} av 

Cly 

by 

kna 

Fiducial 

3.0 

0.648 ± 0.042 

13.1 ±0.42 

1.01 ±0.10 

0.627 ±0.033 

1.66 ±0.013 

4.09 ± 0.040 

2.8 

0.644 ± 0.037 

13.4 ±0.46 

0.979 ± 0.091 

0.610 ±0.030 

1.64 ±0.013 

3.89 ±0.033 

2.6 

0.652 ± 0.034 

13.6 ±0.47 

0.970 ± 0.084 

0.590 ±0.028 

1.61 ±0.011 

3.65 ±0.036 

2.4 

0.666 ± 0.030 

13.5 ± 0.48 

0.963 ± 0.076 

0.561 ±0.027 

1.58 ±0.010 

3.39 ±0.030 

2.2 

0.677 ±0.027 

13.3 ±0.53 

0.961 ±0.07 

0.533 ±0.025 

1.54 ±0.0097 

3.11 ±0.028 

Eulerian 

2.9 

0.719 ± 0.064 

16.1 ±0.73 

0.734 ±0.10 

0.413 ± 0.043 

1.65 ± 0.038 

5.09 ±0.071 

2.6 

0.71 ±0.054 

18.0 ± 1.2 

0.702 ± 0.085 

0.394 ± 0.037 

1.60 ±0.031 

4.76 ± 0.064 

2.3 

0.73 ± 0.046 

19.2 ±1.8 

0.711 ±0.082 

0.370 ± 0.034 

1.55 ±0.028 

4.22 ±0.065 

Lagrange 

3.0 

0.847 ±0.070 

16.4 ±0.95 

1.39 ±0.20 

0.664 ± 0.049 

1.63 ± 0.035 

4.87 ±0.07 

2.8 

0.822 ± 0.063 

17.4 ± 1.3 

1.31 ±0.18 

0.649 ± 0.046 

1.61 ±0.032 

4.70 ± 0.059 

2.6 

0.824 ±0.056 

18.1 ± 1.4 

1.24 ±0.15 

0.619 ±0.043 

1.59 ±0.028 

4.46 ± 0.054 

2.4 

0.836 ± 0.050 

18.1 ± 1.5 

1.16 ±0.13 

0.574 ± 0.038 

1.56 ±0.026 

4.16 ±0.051 

2.2 

0.907 ±0.043 

15.8 ±0.8 

0.981 ±0.085 

0.463 ± 0.029 

1.49 ±0.019 

3.57 ±0.034 

Planck 

3.0 

0.792 ± 0.055 

17.1 ± 1.1 

1.16 ±0.15 

0.578 ± 0.045 

1.63 ± 0.035 

4.79 ± 0.054 

2.8 

0.773 ± 0.050 

19.2 ±1.4 

1.16 ±0.14 

0.608 ± 0.040 

1.65 ±0.031 

4.50 ±0.053 

2.6 

0.781 ± 0.044 

21.1 ± 1.8 

1.15 ±0.12 

0.611 ±0.035 

1.64 ±0.027 

4.19 ±0.051 

2.4 

0.851 ± 0.040 

19.5 ±1.9 

1.06 ±0.093 

0.548 ± 0.029 

1.61 ±0.023 

3.63 ±0.031 

2.2 

0.867 ±0.035 

19.4 ± 1.5 

1.06 ±0.081 

0.514 ±0.027 

1.60 ±0.021 

3.35 ±0.029 
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Table 8. Values of the bias parameters for different physical properties and for 52 = 0. 


z 

P 

brS 

b r7 j 

Fiducial (a = 0.88, log T = 4.3) 

3.0 

1.205 ± 0.049 

0.5546 ± 0.0086 

0.681 ±0.028 

2.8 

1.284 ± 0.052 

0.5577 ± 0.0083 

0.733 ± 0.030 

2.6 

1.343 ± 0.055 

0.5588 ± 0.0084 

0.771 ± 0.032 

2.4 

1.385 ± 0.056 

0.5574 ± 0.0080 

0.796 ± 0.032 

2.2 

1.405 ± 0.061 

0.5536 ± 0.0079 

0.807 ± 0.034 

SO.76 (a = 0.76) 

3.0 

1.086 ± 0.042 

0.6319 ±0.0096 

0.700 ± 0.028 

2.8 

1.161 ±0.045 

0.6373 ± 0.0096 

0.757 ±0.030 

2.6 

1.219 ± 0.047 

0.6409 ± 0.0098 

0.803 ± 0.032 

2.4 

1.257 ±0.050 

0.6428 ± 0.0095 

0.833 ± 0.033 

2.2 

1.284 ± 0.050 

0.6401 ± 0.0094 

0.853 ± 0.033 

SO.64 (a = 0.64) 

3.0 

0.965 ± 0.033 

0.7287 ±0.011 

0.717 ±0.025 

2.8 

1.035 ± 0.035 

0.7366 ±0.011 

0.780 ± 0.028 

2.6 

1.094 ± 0.038 

0.7429 ±0.011 

0.835 ± 0.030 

2.4 

1.137 ±0.039 

0.7478 ±0.011 

0.877 ±0.031 

2.2 

1.163 ±0.041 

0.7490 ±0.011 

0.904 ± 0.032 

G1.3 (7 = 1.3) 

3.0 

1.123 ±0.049 

0.5770 ± 0.0088 

0.661 ±0.029 

2.8 

1.216 ±0.051 

0.5790 ± 0.0089 

0.720 ±0.031 

2.6 

1.290 ± 0.054 

0.5788 ± 0.0086 

0.767 ±0.032 

2.4 

1.348 ± 0.056 

0.5764 ± 0.0079 

0.801 ± 0.033 

2.2 

1.385 ± 0.060 

0.5705 ± 0.0081 

0.820 ± 0.035 

G1.0 (7 = 1.0) 

3.0 

1.090 ± 0.050 

0.5925 ± 0.0091 

0.658 ± 0.030 

2.8 

1.191 ±0.051 

0.5930 ± 0.0086 

0.723 ±0.031 

2.6 

1.275 ± 0.055 

0.5912 ± 0.0084 

0.774 ± 0.034 

2.4 

1.343 ± 0.059 

0.5869 ± 0.0082 

0.813 ±0.035 

2.2 

1.391 ± 0.060 

0.5794 ± 0.0080 

0.836 ± 0.035 

G1T4 (7 = 1.0, log T = 4.0) 

3.0 

1.168 ±0.069 

0.5682 ± 0.0086 

0.677 ± 0.040 

2.8 

1.300 ± 0.071 

0.5679 ± 0.0083 

0.756 ± 0.041 

2.6 

1.425 ± 0.074 

0.5645 ± 0.0084 

0.826 ± 0.042 

2.4 

1.543 ± 0.055 

0.5578 ± 0.0050 

0.887 ±0.031 

2.2 

1.643 ± 0.081 

0.5475 ± 0.0068 

0.934 ± 0.044 
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Table 9. Values of the fitting parameters for different physical properties and for q 2 = 0. 


z 

Qi 

k p [ h / Mpc } 

k [ h / Mpc] av 

a v 

by 

fond 

Fiducial (a = 0.88, log T = 4.3) 


0.648 ± 0.042 

13.1 ±0.4 

1.007 ±0.103 

0.627 ±0.033 

1.66 ± 0.01 

4.09 ± 0.04 

2.8 

0.644 ± 0.037 

13.4 ±0.5 

0.979 ±0.091 

0.610 ± 0.030 

1.64 ± 0.01 

3.89 ± 0.03 

2.6 

0.652 ± 0.034 

13.6 ±0.5 

0.970 ± 0.084 

0.590 ±0.028 

1.61 ± 0.01 

3.65 ±0.04 

2.4 

0.666 ± 0.030 

13.5 ±0.5 

0.963 ± 0.076 

0.561 ± 0.027 

1.58 ±0.01 

3.39 ±0.03 

2.2 

0.677 ±0.027 

13.3 ±0.5 

0.961 ±0.070 

0.533 ±0.025 

1.54 ± 0.01 

3.11 ±0.03 

G1.3 (7 = 1.3) 


0.612 ±0.041 

13.6 ±0.5 

1.044 ±0.115 

0.633 ± 0.036 

1.71 ± 0.02 

4.22 ± 0.04 

2.8 

0.603 ± 0.038 

14.2 ±0.5 

0.996 ±0.103 

0.612 ± 0.034 

1.70 ±0.01 

4.04 ± 0.04 

2.6 

0.609 ± 0.033 

14.6 ±0.6 

0.971 ±0.089 

0.586 ±0.031 

1.67 ±0.01 

3.81 ± 0.03 

2.4 

0.625 ± 0.029 

14.8 ±0.6 

0.951 ±0.080 

0.551 ± 0.029 

1.64 ± 0.01 

3.53 ±0.03 

2.2 

0.641 ± 0.026 

14.7 ±0.7 

0.937 ±0.071 

0.514 ± 0.026 

1.60 ±0.01 

3.23 ±0.01 

G1.0 (7 = 1.0) 


0.611 ±0.042 

13.6 ±0.5 

1.034 ±0.114 

0.614 ±0.036 

1.72 ±0.02 

4.21 ± 0.04 

2.8 

0.600 ± 0.037 

14.2 ±0.5 

0.980 ±0.100 

0.593 ± 0.034 

1.71 ± 0.02 

4.04 ± 0.04 

2.6 

0.605 ± 0.033 

14.7 ±0.6 

0.949 ± 0.090 

0.565 ± 0.032 

1.68 ±0.01 

3.82 ± 0.03 

2.4 

0.621 ±0.029 

14.9 ±0.6 

0.925 ± 0.079 

0.528 ± 0.030 

1.65 ±0.01 

3.55 ±0.02 

2.2 

0.638 ± 0.027 

14.9 ±0.7 

0.909 ±0.071 

0.490 ±0.027 

1.61 ± 0.01 

3.25 ± 0.03 

G1T4 (7 = 1.0, log T 

= 4.0) 






0.620 ± 0.043 

14.4 ±0.5 

0.711 ±0.087 

0.318 ±0.039 

1.60 ±0.02 

4.66 ± 0.06 

2.8 

0.594 ±0.039 

15.5 ±0.7 

0.683 ± 0.079 

0.317 ±0.037 

1.60 ±0.02 

4.53 ±0.06 

2.6 

0.591 ±0.035 

16.5 ±0.9 

0.669 ± 0.070 

0.306 ± 0.034 

1.59 ±0.02 

4.30 ± 0.04 

2.4 

0.604 ± 0.020 

17.1 ±0.9 

0.658 ± 0.044 

0.282 ± 0.020 

1.57 ±0.02 

3.98 ±0.04 

2.2 

0.621 ±0.025 

17.4 ±1.1 

0.651 ±0.058 

0.254 ± 0.028 

1.54 ± 0.02 

3.61 ±0.02 

SO. 76 (a = 0.76) 


0.762 ± 0.055 

13.7 ±0.6 

1.383 ±0.173 

0.788 ±0.041 

1.66 ±0.03 

4.58 ±0.05 

2.8 

0.758 ± 0.050 

14.1 ±0.7 

1.314 ±0.139 

0.763 ± 0.035 

1.64 ±0.03 

4.37 ± 0.04 

2.6 

0.767 ± 0.045 

14.4 ±0.6 

1.270 ±0.122 

0.734 ±0.032 

1.61 ± 0.02 

4.12 ± 0.04 

2.4 

0.783 ± 0.040 

14.6 ±0.7 

1.242 ±0.109 

0.700 ±0.031 

1.58 ±0.02 

3.83 ±0.04 

2.2 

0.804 ± 0.035 

14.5 ±0.6 

1.211 ±0.096 

0.659 ±0.028 

1.55 ±0.02 

3.52 ± 0.03 

SO. 64 (a = 0.64) 


0.958 ± 0.074 

13.6 ±0.6 

2.191 ±0.277 

0.994 ± 0.044 

1.63 ±0.03 

5.17 ±0.05 

2.8 

0.949 ± 0.069 

14.2 ±0.5 

2.015 ±0.238 

0.958 ±0.041 

1.61 ± 0.02 

4.93 ± 0.04 

2.6 

0.958 ±0.061 

14.7 ±0.6 

1.879 ±0.203 

0.915 ±0.038 

1.59 ±0.03 

4.67 ± 0.04 

2.4 

0.979 ± 0.054 

15.1 ±0.8 

1.772 ±0.172 

0.868 ±0.035 

1.56 ±0.02 

4.37 ± 0.04 

2.2 

1.004 ±0.051 

15.3 ±0.8 

1.677 ±0.151 

0.818 ±0.032 

1.53 ±0.02 

4.03 ± 0.04 
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Table 10. Bias parameters when varying F and z independently. 


z, F 

P 

^tS 

b T r) 

Fiducial 

3.0 

1.205 ±0.049 

0.5546 ± 0.0086 

0.681 ± 0.028 

2.8 

1.284 ±0.052 

0.5577 ± 0.0083 

0.733 ± 0.030 

2.6 

1.343 ±0.055 

0.5588 ± 0.0084 

0.771 ± 0.032 

2.4 

1.385 ±0.056 

0.5574 ± 0.0080 

0.796 ± 0.032 

2.2 

1.405 ±0.061 

0.5536 ± 0.0079 

0.807 ± 0.034 

F fixed to 0.781 

3.0 

1.274 ± 0.048 

0.6038 ± 0.0061 

0.790 ± 0.021 

2.8 

1.310 ±0.052 

0.5813 ± 0.0068 

0.782 ± 0.026 

2.6 

1.344 ± 0.056 

0.5588 ± 0.0077 

0.721 ± 0.031 

2.4 

1.383 ± 0.060 

0.5356 ± 0.0093 

0.768 ± 0.039 

2.2 

1.429 ±0.065 

0.5121 ±0.0110 

0.751 ± 0.052 

z fixed to 2.6 

0.696 

1.275 ±0.057 

0.517 ±0.008 

0.676 ± 0.030 

0.740 

1.318 ±0.058 

0.536 ± 0.008 

0.640 ± 0.036 

0.781 

1.343 ±0.057 

0.558 ± 0.008 

0.769 ± 0.032 

0.818 

1.347 ±0.054 

0.581 ± 0.008 

0.804 ± 0.032 

0.852 

1.332 ±0.069 

0.609 ± 0.008 

0.834 ± 0.032 
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Table 11. Non-linear fit parameters when varying F and z independently. 


z 

q\ 

k p [h/Mpc] 

k% v [h/Mpc] av 

Qjy 

by 

^na 

Fiducial 

3.0 

0.648 ± 0.042 

13.1 ±0.4 

1.007 ±0.103 

0.627 ±0.033 

1.66 ±0.01 

4.09 ± 0.04 

2.8 

0.644 ± 0.037 

13.4 ±0.5 

0.979 ± 0.091 

0.610 ±0.030 

1.64 ±0.01 

3.89 ± 0.03 

2.6 

0.652 ± 0.034 

13.6 ±0.5 

0.970 ± 0.084 

0.590 ± 0.028 

1.61 ±0.01 

3.65 ± 0.04 

2.4 

0.666 ± 0.030 

13.5 ±0.5 

0.963 ± 0.076 

0.561 ±0.027 

1.58 ±0.01 

3.39 ±0.03 

2.2 

0.677 ±0.027 

13.3 ±0.5 

0.961 ± 0.070 

0.533 ± 0.025 

1.54 ±0.01 

3.11 ±0.03 

F fixed to 0.781 

3.0 

0.710 ±0.041 

16.0 ±0.8 

1.195 ±0.114 

0.703 ± 0.033 

1.66 ±0.02 

4.12 ±0.04 

2.8 

0.683 ± 0.036 

14.7 ±0.7 

1.070 ± 0.093 

0.644 ± 0.029 

1.64 ±0.02 

3.88 ± 0.04 

2.6 

0.652 ± 0.033 

13.6 ±0.5 

0.970 ± 0.083 

0.590 ± 0.029 

1.61 ±0.02 

3.65 ±0.03 

2.4 

0.620 ±0.030 

12.5 ±0.5 

0.876 ± 0.071 

0.530 ± 0.027 

1.58 ±0.02 

3.43 ± 0.03 

2.2 

0.585 ± 0.028 

11.6 ±0.4 

0.793 ± 0.064 

0.473 ± 0.026 

1.56 ±0.02 

3.22 ±0.04 

z fixed to 2.6 

0.696 

0.585 ± 0.034 

11.7 ±0.3 

0.811 ±0.074 

0.517 ±0.030 

1.62 ±0.02 

3.61 ±0.04 

0.740 

0.614 ±0.034 

12.5 ±0.4 

0.885 ± 0.080 

0.555 ± 0.029 

1.61 ±0.02 

3.65 ± 0.04 

0.781 

0.652 ± 0.034 

13.6 ±0.5 

0.969 ± 0.085 

0.590 ± 0.029 

1.61 ±0.02 

3.65 ±0.03 

0.818 

0.698 ± 0.034 

14.8 ±0.7 

1.063 ± 0.090 

0.620 ± 0.029 

1.61 ±0.02 

3.60 ±0.03 

0.852 

0.751 ± 0.044 

16.4 ± 1.1 

1.160 ±0.118 

0.646 ± 0.035 

1.60 ±0.02 

3.52 ±0.03 











